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In Part I of this study (Ireland et al. 2015), we analyzed the motion of inertial particles in 
isotropic turbulence in the absence of gravity using direct numerical simulation (DNS). 
Here, in Part II, we introduce gravity and study its effect on single-particle and particle- 
pair dynamics over a wide range of flow Reynolds numbers, Froude numbers, and particle 
Stokes numbers. The overall goal of this study is to explore the mechanisms affecting 
particle collisions, and to thereby improve our understanding of droplet interactions in 
atmospheric clouds. We find that the dynamics of heavy particles falling under gravity can 
be artificially influenced by the finite domain size and the periodic boundary conditions, 
and we therefore perform our simulations on larger domains to reduce these effects. We 
first study single-particle statistics which influence the relative positions and velocities of 
inertial particles. We see that gravity causes particles to sample the flow more uniformly 
and reduces the time particles can spend interacting with the underlying turbulence. We 
also find that gravity tends to increase inertial particle accelerations, and we introduce 
a model to explain that effect. 

We then analyze the particle relative velocities and radial distribution functions (RDFs), 
which are generally seen to be independent of Reynolds number for low and moderate 
Kolmogorov-scale Stokes numbers St. We see that gravity causes particle relative ve¬ 
locities to decrease by reducing the degree of preferential sampling and the importance 
of path-history interactions, and that the relative velocities have higher scaling expo¬ 
nents with gravity. We observe that gravity has a non-trivial effect on clustering, acting 
to decrease clustering at low St and to increase clustering at high St. By considering 
the effect of gravity on the clustering mechanisms described in the theory of Zaichik & 
Alipchenkov (2009), we provide an explanation for this non-trivial effect of gravity. We 
also show that when the effects of gravity are accounted for in the theory of Zaichik & 
Alipchenkov (2009), the results compare favorably with DNS. The relative velocities and 
RDFs exhibit considerable anisotropy at small separations, and this anisotropy is quan¬ 
tified using spherical harmonic functions. We use the relative velocities and the RDFs 
to compute the particle collision kernels, and find that the collision kernel remains as it 
was for the case without gravity, namely nearly independent of Reynolds number for low 
and moderate St. We conclude by discussing practical implications of the results for the 
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cloud physics and turbulence communities and by suggesting possible avenues for future 
research. 


1. Introduction 

This is the second part of a two-part paper in which we consider the Reynolds-nunrber 
dependence of inertial particle statistics using direct numerical simulations (DNS) of 
homogeneous, isotropic turbulence (HIT). In Part I of this study (Ireland et al. 2015), 
we used high-Reynolds-number DNS to explore the motion of inertial particles in the 
absence of gravity. We saw that particles with weak inertia preferentially sampled certain 
regions of the turbulence (a phenomenon known as ‘preferential sampling,’ see Maxey 
1987; Squires & Eaton 1991; Eaton & Fessler 1994). By exploring the specific regions of 
the flow contributing to this preferential sampling and using the theory in Chun et al. 
(2005), we were able to understand and model the resulting trends. Particles with higher 
inertia had a modulated response to the underlying turbulence (a phenomenon known 
as ‘inertial filtering’), decreasing the particle kinetic energies and accelerations, and we 
found our DNS data for these quantities to be in excellent agreement the models of 
Abrahamson (1975) and Zaichik & Alipchenkov (2008). Such particles also exhibited 
increased relative velocities and ‘caustics’ (Wilkinson & Mehlig 2005; Wilkinson et al. 
2006), which occur as a result of the particles’ memory of their path-history interactions 
with the turbulence. 

A primary goal of our analysis in Part I was to determine the effect of Reynolds 
number on particle collision rates. It is well-known that droplet growth and precipitation 
in warm, cumulus clouds occurs faster than current microphysical models can predict, 
and the discrepancies are generally linked to turbulence effects (see Shaw 2003; Devenish 
et al. 2012; Grabowski & Wang 2013). We explored droplet motions in turbulence at 
the highest Reynolds numbers to date, and used the results to extrapolate to Reynolds 
numbers representative of those in atmospheric clouds. A secondary motivation was to 
understand the extent to which protoplanetary nebula formation, which depends on 
the collision and coalescence of small dust grains, is affected by turbulence. (A more 
complete explanation of the physical processes involved in cloud and protoplanetary 
nebula formation is provided in Part I.) 

To determine the collision rates, we computed particle radial distributions and relative 
velocities and used the theory of Sundaram & Collins (1997) to calculate the kinematic 
collision kernel from these quantities. We observed that the collision rates of weakly iner¬ 
tial particles (such as those that would be present in the early stage of cloud formation) 
are almost entirely insensitive to the flow Reynolds number. This suggests that particle 
collisions are determined by the small-scale turbulence, and that DNS at low Reynolds 
numbers should capture the essential physics responsible for particle collisions in highly 
turbulent clouds. 

One major simplification in Part I was the neglect of gravitational forces on the par¬ 
ticles. However, as noted in Part I, gravity is not negligible for many particle-laden 
environmental flows. For example, in warm cumulus clouds, the gravitational settling 
speeds of droplets may be an order of magnitude larger than the Kolmogorov velocity, 
suggesting that such droplets fall quickly through the turbulence and may therefore have 
a substantially modified response to the underlying flow (Ayala et al. 2008). Therefore, 
in Part II of our study, we systematically explore the effect of gravity on inertial particle 
statistics. 

The inclusion of gravity, while superficially trivial, adds richness and complexity to 
both the physical interpretation of the results as well as the numerics. Physically, the 
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dimensionality of the parameter space is augmented by one to quantify the strength of 
gravity relative to the turbulence parameters. Here, we introduce a Froude number, Fr, 
which is defined as the ratio of the r.rn.s. of the fluid acceleration to the gravitational 
constant (e.g., see Bee et al. 2014). Even for terrestrial clouds (for which gravity may 
be regarded as a constant), the turbulent acceleration may vary by orders of magnitude 
(see Pruppacher & Klett 1997), leading to significant variations in Fr. To quantify the 
role of gravity, therefore, we must analyze particle statistics over a wide range of St, R\, 
and Fr. In addition to the expanded parameter space, the gravitational vector causes a 
reduction in symmetry for the particle field, from isotropic to axisynnnetric, with the axis 
aligned with the gravitational vector (Ayala et al. 2008; Woittiez et al. 2009; Bee et al. 
2014; Gustavsson et al. 2014). We will investigate the degree of anisotropy as a function 
of the system parameters. Finally, there is a numerical concern that relates to the settling 
speed of the particles. With periodic boundary conditions, particles exiting the bottom 
of the box are reintroduced at the top. If the time required to traverse the box becomes 
smaller than the correlation time of the turbulence, the reflected particles’ preexisting 
correlation with the flow will introduce an unphysical effect (Woittiez et al. 2009). This 
effect can be mitigated by increasing the domain size. We will explore different domain 
sizes to quantify this effect (see Appendix A for a complete discussion). 

Previous DNS of inertial particles subjected to gravity have primarily focused on how 
turbulence alters the terminal velocity (Wang & Maxey 1993; Yang & Lei 1998; Yang & 
Shy 2005; Ireland & Collins 2012; Bee et al. 2014; Good et al. 2014) and the collision 
frequency (Franklin et al. 2007; Ayala et al. 2008; Woittiez et al. 2009; Onishi et al. 
2009; Rosa et al. 2013; Bee et al. 2014). Our work extends those studies by perform¬ 
ing simulations on larger domains, over a wider range of Reynolds numbers, and using 
more particle classes. We also consider the effect of gravity on additional particle statis¬ 
tics (such as fluid velocity gradients, Lagrangian timescales, and particle accelerations) 
to provide new insight into the physical mechanisms responsible for particle collisions, 
and we specifically address the influence of anisotropy on these and other statistics. To 
understand the trends in many of these statistics, we introduce theoretical models and 
compare them with the DNS data. 

The organization of this paper is as follows. In §2, we discuss the numerical methods 
and parameters for our simulations. Single-particle statistics are presented in §3. Within 
this section, we analyze particle velocity gradients (§3.1), mean settling speeds (§3.2), and 
accelerations (§3.3). These statistics provide insight into the mechanisms contributing to 
the clustering and near-contact motion of heavy, inertial particles. To further understand 
these and other relevant collision mechanisms, we directly compute particle-pair statistics 
in §4. We consider the particle relative velocities in §4.1 and radial distribution functions 
in §4.2, and use these data to compute the particle kinematic collision kernel in §4.3. We 
conclude in §5 by summarizing our findings and suggesting some practical implications 
for the cloud physics and turbulence communities. 


2. Overview of simulations 

2.1. Fluid phase 

As in Part I, we perform pseudospectral DNS of HIT on a cubic, tri-periodic domain of 
length £ with N 3 grid points, solving the continuity and Navier-Stokes equations for an 
incompressible flow 


V • u = 0, 


(2.1) 
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Table 1. Simulation parameters for the DNS study. All dimensional parameters are in 
arbitrary units, and all quantities are defined in the text in §2.1 and §2.2. 


Simulation 

I 

II 

III 

IV 

V 

Illb 

Rx 

90 

147 

230 

398 

597 

227 

£ 

167T 

87T 

47T 

2t r 

2tt 

2n 

IS 

0.005 

0.002 

0.0008289 

0.0003 

0.00013 

0.0008289 

e 

0.257 

0.244 

0.239 

0.223 

0.228 

0.246 

i 

1.47 

1.44 

1.49 

1.45 

1.43 

1.43 

i/p 

v! 

55.6 

107 

213 

436 

812 

206 

0.912 

0.914 

0.914 

0.915 

0.915 

0.915 

u'/u v 

4.82 

6.15 

7.70 

10.1 

12.4 

7.65 

^max^7 

1.61 

1.63 

1.68 

1.60 

1.70 

1.67 

N 

1024 

1024 

1024 

1024 

2048 

512 


+ u) x u + V = zA7 2 « + /, (2.2) 

where u is the fluid velocity, u; is the vorticity, p is the pressure, pj is the fluid density, v 
is the kinematic viscosity, and / is a large-scale forcing term that is added to make the 
flow field statistically stationary. In these simulations, deterministic forcing is applied to 
wavenumbers with magnitude k = \[2. Note that the gravity term in the Navier-Stokes 
equation is not shown because it is precisely canceled by the mean pressure gradient, 
and so it has no dynamical consequence on the turbulent flow field. More details of the 
numerical methods can be found in Ireland et al. (2013). 

In Part I, £ = 2n for all the simulations performed. To reduce artificial periodicity 
effects in these simulations with gravity, the domain lengths are here extended to £ = 167r 
(for R x = 90), £ = 8tt (for R x = 147), and £ = 4tt (for R x = 230). ( R x = 2fc v /5/ (Sue) 
denotes the Taylor-scale Reynolds number, where k is the kinetic energy and e is the 
turbulent energy dissipation rate.) The grid spacing is kept constant as the domain size 
is increased, and so the small-scale resolution K max ?7 is approximately constant between 
the different domain sizes (where K max is the maximum resolved wavenumber and rj = 
(y 3 /e) 1 / 4 is the Kolmogorov length scale). In increasing the domain size, we also keep 
the viscosity and forcing parameters the same, and thus both small-scale and large-scale 
flow parameters are held approximately constant. At the two highest Reynolds numbers, 
we expect periodicity effects to be minimal, and the domain sizes are the same (£ = 2n) 
both with and without gravity. Refer to Appendix A for a detailed examination of the 
effect of the domain size on fluid and particle statistics. 

The simulation parameters are given in table 1. The simulation results from this study 
will be frequently compared to those from Part I. These comparisons are justified, since 
in all cases, the parameters are very close to those used in Part I, and based on the 
results in Appendix A, we can safely assume that the differences in the particle statistics 
between these simulations and those in Part I are due entirely to gravitational effects and 
not to any differences in the underlying flow fields. In showing comparisons between the 
present data and those from Part I, we will refer to both fields by the Reynolds numbers 
from table 1, keeping in mind that the corresponding Reynolds numbers from Part I may 
differ slightly (but by no more than 5%). 

To perform a more complete parametric study of the effects of inertia and gravity 
on particle statistics, we also conducted a simulation with similar flow parameters to 
simulation III, but with a smaller domain size (due to computational limitations). The 
parameters for this simulation (referred to as Illb) are also given in table 1. 
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2.2. Particle phase 

We simulate the motion of spherical particles with finite inertia and gravitational forces. 
To model the dynamics of inertial particles, we make the following simplifying assump¬ 
tions. The particles are assumed to be small (d /77 U where d is the particle diam¬ 
eter) and dense ( p p /pf 1, where p p is the particle density), and subject to only 

linear drag forces. The last assumption is reasonable when the particle Reynolds num¬ 
ber Re p = | u p (t) — v p {t)\/v < 0.5, where u p (t) = u(x p (t),t) denotes the undisturbed 
fluid velocity at the particle position x p , and v p denotes the velocity of the particle (El- 
ghobashi & Truesdell 1992). (As in Part I, we use the superscript p on x, u, and v to 
denote time-dependent, Lagrangian variables defined along particle trajectories. Phase- 
space positions and velocities will be denoted without the superscript p.) In a recent 
study, Good et al. (2014) found that linear-drag models yield settling speeds that are 
inconsistent with experiments at larger settling velocities, demonstrating the breakdown 
of the linear-drag model for St > 3. However, simply introducing a nonlinear drag coeffi¬ 
cient does not accurately capture the nonlinear drag effects in a time-dependent flow like 
turbulence. A full wake-resolving, nonlinear model would be required, but is beyond the 
scope of this study. Despite the need for caution in the interpretation of the higher-Stokes 
simulations, linear drag yields the correct qualitative trends, thus making the physical 
arguments and discussions presented in this paper valid even if a more accurate nonlinear 
model were introduced. Moreover, the added value of using linear drag is that the major¬ 
ity of theoretical models make the same assumption, and so this facilitates comparisons 
between DNS and theory. 

Under these assumptions, the governing equations for the inertial particles are (Maxey 
& Riley 1983) 

d 2 x p dv p u(x p (t),t) - v p (t) 

7ir = lu = - 7„ -+ »• (23) 

where g = (0,0,— g) is the gravitational acceleration vector and r p = p p d 2 /(l&pfv) is 
the response time of the particle. To compute u(x p (t),t), we employ an eight-point B- 
spline interpolation from the Eulerian grid (van Hinsberg et al. 2012). As in Part I, we 
begin computing particle statistics once the particle distributions and velocities become 
statistically stationary and independent of their initial conditions. For a subset A/ rac k e d 
of the total number of particles in each class N p , we store particle positions, velocities, 
and velocity gradients every O.lr^ for a duration of about 100r^. These data are used to 
compute Lagrangian correlations, accelerations, and timescales. 

The Stokes number St = t p /t v is a non-dimensional measure of the particle inertia, 
where t p = (v/e) 1 / 2 is the Kolmogorov timescale. Gravity can be parameterized in a 
number of ways (e.g., see Good et al. 2014). Here, we define the Froude number as 
Fr = e 3 / 4 /(iz 1 / 4 g), where e 3 / 4 /iA/ 4 is the Kolmogorov estimate for the acceleration r.m.s. 
(Note that the Froude number is sometimes expressed as the inverse of our definition; we 
choose the present definition for consistency with the standard practice of defining the 
Froude number with g in the denominator.) Alternatively, gravity can be parameterized 
by the settling parameter Sv = t p q/u v , which is the ratio of the particle’s settling velocity 
T p g to the Kolmogorov velocity of the turbulence u v = {ve) 1 ^. Note that these two non- 
dimensional parameters are related through the Stokes number, Fr = St/Sv, and that 
the Froude number Fr is independent of r p . 

Since a primary aim of this paper is to study the effect of gravity at conditions repre¬ 
sentative of those in cumulus clouds, we calculate Fr and Sv by assuming a dissipation 
rate e = 10~ 2 m 2 /s 3 (Shaw 2003), a kinematic viscosity v = 1.5 x 10 -5 m 2 /s, and a 
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gravitational acceleration g = 9.8 m/s 2 , yielding Fr = 0.052 or Sv = 19.3St. Simulations 
I, II, III, and IV were run with Fr = 0.052 and particles with Stokes numbers in the 
range 0 ^ St ^ 3; due to computational limitations, simulation V was only run with 
0 ^ St ^ 0.3. 

However, experimental observations suggest that cloud dissipation rates can vary by 
orders of magnitude (e.g., see Pruppacher & Klett 1997), resulting in commensurate 
variations in Fr and Sv. For example, a strongly turbulent cumulonimbus cloud with 
e ~ 10 -1 m 2 /s 3 yields Fr « 0.3 and Sv ss 3ASt, while a weakly turbulent stratiform 
cloud with e ~ 10 -3 m 2 /s 3 yields Fr ~ 0.01 and Sv ~ lOOSt (Pinsky et al. 2007). 
To study a larger St-Sv parameter space, we analyzed 513 different combinations of St 
and Sv in simulation Illb, with 0 ^ St ^ 56.2 and 0 ^ Sv ^ 100. The results from 
this simulation will be used to study detailed trends in particle accelerations, clustering, 
relative velocities, and collision rates for different values of particle inertia and gravity. 
(The trends in the particle kinetic energies and settling velocities obtained from this data 
set are discussed in detail in Good et al. (2014).) 

3. Single-particle statistics 

As in Part I, we begin with the statistics of individual inertial particles. We present 
velocity gradient statistics in §3.1, gravitational settling velocity statistics in §3.2, and 
acceleration statistics in §3.3. 

3.1. Velocity gradient statistics 

We begin by considering velocity gradients sampled along inertial particle trajectories, 
here denoted by A(x p (t),t) = Vu(x p (t),t). These statistics affect the relative motion of 
particles with weak-to-moderate inertia with separations in the dissipation range, and 
also quantify the degree to which the inertial particles preferentially sample the un¬ 
derlying turbulence (e.g., see Chun et al. 2005). This information will prove useful in 
interpreting the relative velocity and clustering statistics presented in §4.1 and §4.2. It 
is convenient to define the rate-of-strain and rate-of-rotation tensors as S(x p (t),t) = 
[A(x p (t),t) + A T (x p (t),t)\/ 2 and 7 Z(x p (t),t.) = [A(x p (t),t) — A 1 (x p (t) 1 t)\/2, respec¬ 
tively. 

We first consider the variance of the diagonal components of S averaged along the par¬ 
ticle trajectories. For fluid particles that uniformly sample the flow, their variances reduce 
to ([Sn(a: p (t),f)] 2 ) = ([. S 2 2 (x p (t) 1 t )] 2 ) = ([S 33 (x p (t), t)] 2 } = l/(15r 2 ). In isotropic tur¬ 
bulence without gravity, these components are statistically equivalent; however, in the 
presence of gravity, the symmetry of the particle motion is reduced from isotropic to ax- 
isynnnetric, and thus the statistics of S along the direction of gravity, x 3l will differ from 
the transverse directions, x 3 and X 2 - For the remainder of the paper, the x 3 direction 
will be denoted as the ‘vertical’ direction and the x 3 and x 2 directions will be referred to 
as the ‘horizontal’ directions. The data for the two horizontal directions will be averaged 
together whenever possible, and denoted with the subscript ‘1.’ 

Figure 1 shows the variance of (a) horizontal ((Sf 1 ) p = ([<Sn(£e p (f),f)] 2 )) and (b) 
vertical (( S 33 ) p = ([S 33 (x p (t),t)] 2 )) components of the rate-of-strain tensor as a function 
of St at different values of R\, both with gravity (Fr = 0.052) and without gravity 
(Fr = 00 ). Without gravity, the strain rates decrease with increasing St for St < 0.3. 
This is closely related to the trend in the strain rates observed in Part I, which was 
attributed to the fact that inertia causes particles to be ejected from vortex sheets. With 
gravity, the strain rates in the horizontal direction also decrease with increasing St for 
St < 0.3, as seen in figure 1(a). In addition, they are quite close to the corresponding 
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Figure 1. The normalized variance of the longitudinal velocity gradients sampled by inertial 
particles for different values of R\ and St. Open symbols denote data with gravity (Fr = 0.052), 
and filled symbols denote data without gravity (Fr = oo). The horizontal dotted line indicates 
the expected value for uniformly distributed particles in isotropic turbulence. The gradients in 
the horizontal and vertical directions are shown in (a) and (b), respectively. 


values without gravity, suggesting that gravity does not lead to significant changes in 
preferential sampling in the horizontal direction at low St. We see in figure 1(b), however, 
that gravity causes low-St particles to sample regions of larger vertical strain rates, and 
that these strain rates are considerably different from those when gravity is absent. We 
also observe that gravity reduces the degree of preferential sampling at high St, causing 
the strain rates to approach the values predicted for uniformly distributed particles in 
isotropic turbulence for St > 1. 

In figure 2, we plot t 2 (S 2 ) p , t 2 (TZ 2 ) p , and t 2 ((S 2 ) p —(TZ 2 ) p ), where ( S 2 ) p = (S(x p (t),t) : 
S(x p (t),t)) and (7 Z 2 ) p = (lZ(x p (t),t) : lZ(x p (t), t)) are the second invariants of the rate- 
of-strain and rate-of-rotation tensors, respectively, for particles with (Fr = 0.052) and 
without (Fr = oo) gravity. As noted in Part I, when preferential sampling effects are 
absent, t 2 (S 2 ) p = t 2 (IZ 2 ) p = 0.5. We see from figure 2 that gravity reduces the degree 
of preferential sampling, causing t 2 (S 2 ) p , t 2 (It 2 ) p , and t 2 ((S 2 ) p — (If 2 ) p ) to be closer 
to the corresponding values for uniformly distributed particles. We also note that prefer¬ 
ential sampling effects are eliminated altogether for St > 1, which is consistent with the 
results presented earlier. The trends in the mean strain and rotation rates with R\ are 
similar both with and without gravity, and are discussed in Part I. 

To explain the reduction in preferential sampling with gravity, we note that gravity, by 
causing particles to settle relative to the underlying flow, reduces the interaction times 
between the particles and the turbulent eddies. As a result, particles have less time to be 
affected by the straining and rotating regions of the flow, and therefore experience less 
preferential sampling (see also Gustavsson et al. 2014). 

We test this argument by computing Lagrangian rate-of-strain and rate-of-rotation 
timescales along inertial particle trajectories, both with and without gravity. These 
timescales will prove useful for extending the theoretical model of Zaichik & Alipclrenkov 
(2009) to account for gravitational effects. Their model, developed for an isotropic parti¬ 
cle field, uses a single timescale for the rate-of-strain tensor and a single timescale for the 
rate-of-rotation tensor. To generalize the model for the case with gravity, directionally 
dependent times scales are required. We discuss the generalized model and its comparison 
with DNS in §4.2. 
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Figure 2. (S 2 ) p (a), {1Z 2 ) P (b), and (S 2 ) p — {TZ 2 ) P (c) as function of St for different values 
of R\. Open symbols denote data with gravity (Fr = 0.052), and filled symbols denote data 
without gravity (Fr = oo). 


We consider both the Lagrangian strain timescale Tjf , which we define as 


rpp 

JL o. . C — 

^ij ^ km 


(Sij(x P ( 0), 0 )Skm(x p (s), s)) ds 


('Sij (®P(0), 0)5 fcm (x p (0), 0)) 


and the Lagrangian rotation timescale T^. , defined as 


(3.1) 


J (n ij (x p (o),o)n km (x p (s),s))ds 

= (K XJ (x p (0),0)K krn (xP(0),0)) 


(3.2) 


In both of these equations, no summation is implied by the repeated indices. 

In the absence of gravity, the particle field is isotropic, which implies the rate-of-strain 
timescales T% 9 ,T% o , T? <- ,1% 9 , T% o , TS <- , Tt 9 , T? o , and To o 

7 *->11 < - > 22 7 *>11^33 7 *>12‘->12 7 ^lS^lS 7 022‘>22 7 *>22«>33 7 <->23‘- > 23 7 

are statistically equivalent, as are the rotation timescales T^ ioKi2 , ^ , and ^ . 
However, with gravity, the nine rate-of-strain timescales are no longer equivalent, nor are 
the three rate-of-rotation timescales. Applying symmetry analysis, we can group the nine 
rate-of-strain timescales into six categories, and three rate-of-rotation timescales into two 
categories, as shown in table 2. Additionally, we introduce Tg S and T^ K as the average 
of the nine rate-of-strain timescales and three rate-of-rotation timescales, respectively. 
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Table 2. Categorization of the strain and rotation timescales with gravity. (Gravity is aligned 
with the * 3 -direction.) The elements within a given row are statistically equivalent, while the 
elements in different rows may differ due to the anisotropy in the particle motion. 


( 1 ) 

rnP 

1 S 11 S 11 

rnP 

*522*522 

(2) 

rnP 

^ 1 1*522 


(3) 

rnP 

± «5l2*5l2 


(4) 

rrp 

*513*513 

r r P 

*523*523 

(5) 

rnP 

*5i 1*533 

rnP 

*522*533 

(6) 

rnP 

± *533*533 


( 1 ) 

rnP 

- L 'Rl2Tli2 


(2) 

rnP 

^13^13 

rnP 

7^23^23 




Figure 3. The averaged strain timescale T^ s (a) and the averaged rotation timescale (b) 
plotted as a function of St for different values of R\. Open symbols denote data with gravity 
(Fr = 0.052), and filled symbols denote data without gravity (Fr = oo). 


Figure 3 shows the averaged time scales Tg S and as a function of the particle 
Stokes number at different values of R\, with and without gravity. In agreement with the 
physical explanations offered above, we see that the timescales without gravity always 
exceed those with gravity. We also note that both with and without gravity, the strain 
timescales are almost entirely insensitive to the Reynolds number. In addition, while the 
rate-of-rotation timescales without gravity vary weakly with R\ (as noted in Part I), 
they are nearly independent of the Reynolds number with gravity. 

We now consider the anisotropy in the different timescales with gravity at R\ = 398 
in figure 4. In this plot, we have averaged together the statistically equivalent timescales 
grouped in table 2, and have denoted these timescales by the first (or only) item in 
these groups. As expected, the timescales show considerable anisotropy with gravity. If 
we consider the limit of strong gravity (Sv v!/u v ), we can derive a simple model for 
the timescales (see Appendix B for the derivation and estimates of droplet diameters 
required for the model to be valid in different atmospheric clouds). Using the relations 
in §B.l, the Lagrangian rate-of-strain and rate-of-rotation timescales take the form 

(■SijSk m ,3 


rpP 

1 SijSk 


Sv 


(3.3) 
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Figure 4. The Lagrangian strain (a) and rotation (b) timescales in different directions at 
R\ = 398, plotted as a function of St, for particles with gravity (Fr = 0.052). 


and 


T p 

HijUkn 


^IZijlZkn 


(3.4) 


respectively, where we have used the top-hat symbol to denote a variable normalized 
by Kolmogorov scales, and Is i:i S km ,3 (^, 7 ^™, 3 ) is the integral lengthscale of SijS km 
(HijUkm) evaluated along the a^-direction. 

Equations (3.3) and (3.4) imply in the limit Sv v!/u v , the Lagrangian integral 
timescales of the flow experienced by the particles are directly proportional to the Eule- 
rian integral lengthscales of the strain and rotation fields. We can therefore write each of 
the Lagrangian timescales in table 2 as functions of their corresponding Eulerian integral 
lengthscales. This allows us to use tensor invariance theory and the properties of the 
strain and rotation rate tensors in HIT to predict the relative magnitude of the different 
timescales. We therefore can show that 


and 


'T'P _ 'T'P ^ 'T'P ^ 'T'P _ 'T'P / 'T'P 

»5ii«S33 ^ 33^33 ‘Sn«Sn ‘Si2<Si2 ‘Si3«Si3 ‘Sll«S22’ 


rpP / rpp 


(3.5) 

(3.6) 


The DNS data in figure 4 follow the trends predicted by (3.5) and (3.6) at high St. 


3.2. Mean particle settling velocities 

We now analyze the effect of turbulence on the mean settling speed of inertial particles. 
The role of the mean settling velocity, aside from its importance to the process of cloud 
rain out, is also related to the evolution of the droplet size distribution through its effect 
on the collision kernel of different size particles, as discussed in Davila & Hunt (2001) 
and Ghosh et al. (2005), and later modeled by Grabowski & Wang (2013). 

Recently, Good et al. (2014) measured settling speeds of particles in turbulence and 
showed that turbulence causes particles with low-to-moderate Stokes numbers to settle 
more quickly than they would in a quiescent fluid. The so-called ‘preferential sweeping’ 
or ‘fast tracking’ of particles (particles sampling downward-moving fluid over upward- 
moving fluid) had been observed in DNS (Wang & Maxey 1993; Yang & Lei 1998; Ireland 
& Collins 2012; Bee et al. 2014) and earlier experiments (Aliseda et al. 2002; Yang & Shy 
2003, 2005). For St > 1, Good et al. (2014) observed a transition to so-called ‘loitering’ 
(particles settling more slowly than they would in a quiescent fluid), again consistent with 
earlier experiments (Nielsen 1993; Yang & Shy 2003; Kawanisi & Shiozaki 2008), although 
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Figure 5. Turbulence-induced enhancements (positive) or reductions (negative) in the mean 
settling velocities of inertial particles with gravity (Fr = 0.052), normalized by the Kolmogorov 
velocity u v . The symbols denote different values of R\. 

Good et al. (2014) were the first to measure it for droplets in air. In contrast, DNS 
with linear drag shows preferential sweeping at all Stokes numbers, with no transition to 
loitering (Wang & Maxey 1993; Yang & Lei 1998; Ireland & Collins 2012; Bee et al. 2014). 
In order to observe loitering in DNS, the linear drag must be replaced by a nonlinear 
drag model (e.g., Clift et al. 1978). The nonlinear drag model yields settling speeds 
that are only in qualitative, but not quantitative, agreement with the experiments, most 
likely due to the simplistic nature of the nonlinear drag model that was used (Good et al. 
2014). 

Here, we extend the range of Reynolds numbers for 0 ^ St ^ 3 and Fr = 0.052 using 
a linear drag model. From (2.3) we can show that 

(u(x p (t), t)) = (v p {t)} - t p q = -(Av) p , (3.7) 

where t p q = (0,0, — r p g) T is the gravitational settling velocity in a quiescent flow. Pref¬ 
erential sweeping corresponds to (Au 3 ) p > 0 and loitering to (Ar; 3 ) p < 0. Figure 5 shows 
(Ai> 3 ) p , normalized by the Kolmogorov velocity u v , as a function of the particle Stokes 
number. Consistent with the DNS results presented in Good et al. (2014), the DNS with 
linear drag shows no statistically significant loitering at any Reynolds number. Moreover, 
we see that the mean settling speeds are independent of R\ for St f 0.1, suggesting that 
in this limit, they are determined entirely by the small-scale turbulence, in agreement 
with Bee et al. (2014). At higher St, the settling speeds are stronger functions of the 
Reynolds number; however, these results must be interpreted with some caution given 
the concerns raised about linear drag in this regime. Nevertheless, we see at least the po¬ 
tential for a significant Reynolds number effect on the mean settling velocities at higher 
Stokes number. 


3.3. Particle accelerations 

We now consider the acceleration statistics of the particles. As noted in Chun et al. 
(2005), the relative velocity of different-sized particles (and thus their collision rate) is 
related to their accelerations. Figure 6 shows the particle acceleration variances ( a\) p = 
((dv p (t)/dt) 2 ) and (a§) p = ((dv p (t)/dt) 2 ), with and without gravity. Notice that the 
acceleration variances with gravity always exceed those without gravity, in some cases by 
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Figure 6. The variance of the particle accelerations in the horizontal (a) and vertical (b) 
directions. Open symbols denote data with gravity (Fr = 0.052), and filled symbols denote data 
without gravity (Fr = oo). 


as much as an order of magnitude. Similar trends were observed in computational and 
experimental studies of channel flows (Gerashchenko et al. 2008; Lavezzo et al. 2010), 
and more recently in the computational study of Parishani et al. (2015) in HIT. 

We also see that gravity reverses the trends in the particle accelerations at low St, caus¬ 
ing the accelerations to increase with increasing St. This trend was recently reported in 
Parishani et al. (2015) at lower values of R\ (R\ = 119 — 143). Our results in Appendix A, 
however, suggest that the acceleration statistics in Parishani et al. (2015), while quali¬ 
tatively correct, have quantitative errors due the influence of artificial periodicity effects 
on the small domain sizes in their simulations. Experimental studies of inertial particle 
accelerations in HIT (Ayyalasomayajula et al. 2006) and Von Karman flows (Volk et al. 
2008a, b) have only observed a monotonic decrease of the particle accelerations with in¬ 
creasing St, presumably because the Froude numbers in those studies (Fr ~ 0.3 for HIT 
and Fr ~ 30 for the Von Karman flow) were too large to observe this trend reversal. 
(Recall that in our study, Fr = 0.052.) 

We now provide a physical explanation for these large accelerations. We begin by 
considering the limit St <C 1. If gravitational forces are negligible, to leading order the 
particle acceleration a p (t) is equivalent to the fluid acceleration at the particle location 
(e.g., see Bee et al. 2006) 

a p (t) ~ t) _|_ u ( x P(t^ £) . \/u(x p (t), t). (3.8) 

As the gravitational force increases (i.e., as Sv increases), it is the interplay between grav¬ 
ity and turbulence that determines the particle acceleration. For example, the downward 
drift due to the gravitational settling will cause the particle to experience different regions 
of the turbulence. We refer to this phenomenon as the ‘gravitational trajectory effect.’ 
(Note that Yudine (1959) coined ‘crossing trajectories’ to refer to the same mechanism.) 

The particle acceleration in the limit of low St and high Sv (i.e., the limit Fr —> 0) 
can be approximated by the derivative of the fluid velocity along the inertial particle 
trajectory (denoted here as du p /dt), yielding (see Bee et al. 2006) 

a'M » f ^ + v , (t) . Vo(x»(t,,<). (3.9) 

We compare the variances of the horizontal and vertical components of du p (t)/dt (de¬ 
noted as (( du\/dt) 2 ) p and ((du 3 / dt) 2 ) p , respectively) with ( a\) p and (a|) p in figure 7. 
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Figure 7. Inertial particle acceleration variances in the horizontal (a) and vertical (b) directions 
for different values of St and R\. All data are with gravity included (Fr = 0.052). Open symbols 
denote ( a\) p (a) and {af) p (b) and filled symbols denote {(dui/dt) 2 ) p (a) and {(du 3 /dt) ) p (b). 


The variances of a p (t) and du p {t)/dt are almost identical at low St, as expected. As 
St increases, the particles have larger settling velocities, causing the fluid seen by the 
particles to change more rapidly in time. As a result, the variance of du p /dt increases 
monotonically without bound with increasing St. In contrast, the inertial particle ac¬ 
celeration peaks around St ~ 1 and then decreases, causing the variance of du p /dt to 
exceed that of a p (t) for St > 0.3. The trend for the acceleration variance occurs because 
the particle’s inertia also modulates its response to the underlying fluid through ‘inertial 
filtering’ (Bee et al. 2006; Salazar & Collins 2012). The particle is unable to respond 
to fluid fluctuations with frequencies greater than 0(l/r p ), causing a reduction in the 
acceleration variance. Beyond St ~ 1, inertial filtering dominates the enhancement due 
to gravitational settling, causing the particle accelerations to decrease with increasing 
St. 

Parishani et al. (2015) also observed a peak in the particle acceleration variance for 
St ~ 1 with gravity, and argued that the preferential sweeping mechanism identified in 
Wang & Maxey (1993) is responsible for producing these large accelerations. While we 
observe a similar peak in the acceleration variances, we find preferential sweeping to be 
negligible at St ~ 1 for R\ < 230 (see §3.2). Our data also indicate that the preferential 
sweeping of St ~ 1 particles varies significantly with R\, while the acceleration variances 
are nearly independent of R\. These findings suggest that preferential sweeping cannot 
fully explain the trends in the acceleration variances with gravity. In contrast, the ap¬ 
proach we have taken is able to predict the dependence of the accelerations on St and 
Sv without appealing to preferential-sweeping arguments. 

By comparing figure 7(a) and figure 7(b), we notice that ( a\) p generally exceeds ( af) p . 
We use the acceleration model derived in §B.2 to explain this result in the limit Sv 
u'/urf. From §B.2, we have the following models for the particle acceleration variances in 
the high-Su limit: 


v! 1 2 T Sul 

2 St.Sv + 3 

u v \ L-stJ 

2 (stSv+^Y 


(3.10) 


(«j) p _ 

a v l u v. 



and 


2 


Sv 

si 


(3.11) 
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Figure 8. Comparision between DNS data and theoretical predictions for the variance of the 
particle accelerations in the horizontal (a) and vertical (b) directions at high St. Open symbols 
denote DNS data with gravity (Fr = 0.052), and the lines in (a) and (b) indicate the predictions 
from (3.10) and (3.11), respectively, for the different Reynolds numbers simulated. 


By taking the ratio of (3.10) and (3.11), we see that 

<«;>- 2S ‘ s » +3 0i) 

<“!>' 2StSv + 2 (|)' ' ’ 

and thus the model also predicts that ( a\) p exceeds (a|) p by an amount that increases 
with R\. 

In the derivation of this model, we assumed that the particles’ primary motion is in 
the vertical direction, from which it followed that the particle accelerations are related to 
the correlation lengtlrscales of the fluid velocities sampled along the particle trajectories. 
Since in isotropic turbulence the longitudinal integral lengthscale is twice the transverse 
integral lengthscale, the particles (as they fall vertically) will experience vertical fluid 
velocities that are more strongly correlated over a time t p than are the horizontal fluid 
velocities. (This phenomenon was referred to as the ‘continuity effect’ in Csanady (1963), 
and is also explained in Yudine (1959) and Good et al. (2014).) As a result, the horizontal 
fluid velocities sampled by the particles will change more rapidly, leading to larger particle 
accelerations in those directions. 

In Parishani et al. (2015), the authors argued that the differences in the horizontal and 
vertical particle accelerations are linked to the preferential sweeping mechanism of Wang 
& Maxey (1993). That is, a falling particle will accelerate toward the downward-moving 
part of an eddy, where the vertical fluid velocity changes slowly, but the horizontal 
fluid velocity changes rapidly. While this explanation is plausible for St < 0.3 (where 
preferential sweeping is strong, as shown in §3.2), it fails to explain the strong difference 
between the horizontal and vertical accelerations at higher values of St, where we observe 
negligible preferential sweeping for R\ < 230. The models developed by Parishani et al. 
(2015) are for the low -St limit, whereas our approach is based on the characteristics of 
the Eulerian velocity held and provides a quantitative prediction for the accelerations at 
larger values of St. 

In figure 8, we compare the DNS data for ( a\) p and (a§) p to the model predictions for 
St^l (Sv ^ 19.3). We find that the results are in good agreement at the largest values 
of St, with the model being able to predict both the trends with R\ and the decrease in 
the variances with increasing St. 
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Figure 9. Filled contours of the particle acceleration variances, ( ai) p (a) and (a§) p (b), normal¬ 
ized by the Kolmogorov acceleration variance a )), for R\ = 227. The contours are logarithmically 
scaled, and the labels on the colorbar denote exponents of the decade. The diagonal lines de¬ 
note three different values of Fr, corresponding to conditions representative of stratiform clouds 
(Fr = 0.01), cumulus clouds (Fr = 0.05), and cumulonimbus clouds (Fr = 0.3). The dotted 
lines corresponds to (oi) p = a^ in (a) and (a§) p = in (b). 



To further explore the trends in the accelerations for particles with varying levels of 
inertia and gravity, we plot in figure 9 the inertial particle acceleration variances for 
0 < St < 56.2, 0 < Sv < 100, and R\ = 227. (Recall that statistics at the largest 
values of Sv have errors due to unphysical effects of the periodic boundary conditions. 
Refer to Appendix A for a more detailed explanation.) As expected, the particle ac¬ 
celeration variances are the largest for St C 1 and Sv 1, due to the gravitational 
trajectory effect, as predicted by (3.10) and (3.11). For sufficiently large St, particles 
filter out nearly all of the large-scale turbulence, and so acceleration variances are small. 
For intermediate St and Sv, the particle accelerations are determined by a combination 
of preferential-sampling, gravitational-trajectory, and inertial-filtering effects. We also 
observe from figure 9 that the particle acceleration variances are the largest under condi¬ 
tions representative of stratiform clouds (Fr ~ 0.01) and the smallest under conditions 
representative of cumulonimbus clouds (Fr ~ 0.3). 


4. Two-particle statistics 

We now consider particle pair statistics, in particular, particle relative velocities (§4.1), 
clustering (§4.2), and the collision kernels (§4.3). In each case, we compare our results to 
those without gravity (from Part I) to highlight the role of gravity. 

4.1. Particle relative velocities 

In this section, we examine the effect of gravity on the relative velocities of two identical 
inertial particles. We first discuss the expected effect of gravity on the relative velocities 
from a theoretical framework (§4.1.1), and then analyze the DNS results (§4.1.2). 

4.1.1. Theoretical framework for particle relative velocities 

As in Part I, we define the relative position and velocity as r p (t ) and w p (t), respec¬ 
tively, and Ait (r p (t),t) as the difference in fluid velocities at the particle locations. By 
subtracting (2.3) for two particles at a given instant in time, we get the following equation 
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for the relative position and velocity of the two particles 

d 2 r p dw p A u(r p (t),t) — w p (t) 
dt 2 dt t p 

Notice that the gravity vector precisely cancels out of this equation, making it appear as 
though gravity does not influence the relative motion of the particle pair; however, this 
is not true, as gravity has an implicit effect on the statistics of Art (r p (t),t), and this 
substantially modifies the relative position and velocity statistics of the particles (as will 
be shown below). 

Following the nomenclature in Part I, we define the second-order particle relative 
velocity structure function tensor as 

Sf(r) = (w p (t)w p (t)^, 


where (•),. denotes an ensemble average conditioned on r p (t) = r. We now use the 
formal solution of (4.1) to construct the exact expression for S^r). This will allow us 
to explain the expected trends in the structure functions with inertia and gravity. The 
exact expression for Sf (r) is given by (e.g., see Pan & Padoan 2010) 


S P 2 (r) 



Aw(r p (s), s)Au(r p (S), Sn ^ exp St 1 (s + 5) 


dsdS , 


(4.2) 


where Z denotes a variable Z normalized by Kolmogorov scales. 

We first consider (4.2) in the absence of gravity. From (4.2), we observe that particle 
relative velocities are affected by the fluid velocity difference A.u(r p (t),t.) along their 
trajectories. For small values of St, the high damping of the exponential term in (4.2) 
causes contributions to the integral from the particle history to be negligible, and hence 
the particle relative velocities in the equation become equivalent to the instantaneous 
fluid velocity difference at the current particle position, reflecting the fact that preferen¬ 
tial sampling is the dominant mechanism affecting relative particle motion. As shown in 
Part I, preferential sampling causes inertial particles to experience larger (smaller) fluid 
velocity differences than those of fluid particles in directions parallel (perpendicular) to 
the particle separation vector. As St increases, in the absence of gravity, the exponential 
damping diminishes and the influence of the particle history becomes increasingly impor¬ 
tant. At higher St, particles retain a memory of their interactions with the turbulence, 
and it is the fluid velocities at larger separations along their path histories that dominate 
the contribution to the particle velocity dynamics, leading to particle relative velocities 
that are larger than the local fluid velocity difference. This phenomenon has been referred 
to as ‘caustics’ (e.g., see Wilkinson & Mehlig 2005; Wilkinson et al. 2006). Refer to Bragg 
& Collins (20146) for a more complete theoretical description of this phenomenon, and 
to Part I for supporting DNS evidence. 

We now introduce gravity and discuss its effect on the relative velocities. Recall in 
§3.1 that gravity generally reduces the degree of preferential sampling. With gravity, 
we therefore expect the relative velocities of small-S 1 ! particles to be affected less by 
preferential sampling, and therefore to have dynamics closer to those of fluid particles. 
At higher values of St, gravity alters the inertial particle relative velocities through its 
influence on the path-history mechanism. We observed in §3.1 that gravity reduces the 
correlation timescales of the flow along particle trajectories. This implies that gravity 
causes the fluid velocity differences in (4.2) to become decorrelated more rapidly, and 
thus it reduces the correlation radius over which the relative velocities are influenced by 
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their path-history interactions with the fluid. We therefore expect the particle relative 
velocities to be reduced by gravity. 

In addition to the magnitude changes caused by gravity, its directional nature causes 
the symmetry of the particle field to become anisotropic. We define the mean inward, 
parallel and perpendicular relative velocity statistics as 



,° 

/ w\\p(w\\\r)dw\\ , 

J —OO 

(4.3) 


/ H] P{w\\\r)dw^ , 

/ — OO 

(4.4) 

W=- 

f°° 2 

/ Kl] p(w±\r)dw± , 

/ — OO 

(4.5) 


where and w± are the parallel and perpendicular components of the relative velocity, 
respectively, and p(w^\r) and p{w±\r) are their respective probability density functions, 
conditioned on the separation vector r. In the absence of gravity, Sff , S^l, and S'^y are 
functions of the separation distance r = \r\ only. With gravity, these statistics become 
functions of the angle between the separation vector r and the gravity vector g, here 
denoted as 6. A convenient way to express this dependence is to expand these functions 
in a spherical harmonic series 


OO 


St ||W = E C '2°t(r)^W , 

(4.6) 

e- o 


^||(r) = f;4W^(«) , 

(4.7) 

t=0 


SlAr) = T,E° u {r)Y^{e) , 

t=o 

(4.8) 


where C^r), D^r), and E^fr), are the spherical harmonic coefficients, and F 2 ° are the 
spherical harmonic functions. In discussing the results below, we will quantify the degree 
of anisotropy introduced by gravity by analyzing the relative magnitudes of the spherical 
harmonic coefficients. Similar expansions will be used to describe the anisotropic radial 
distribution function that results from gravity. 

4.1.2. Relative velocity results 

We now analyze the results from the DNS using the decompositions described in §4.1.1. 
To study how gravity affects the magnitude of the relative velocity, we first consider the 
zeroth-order coefficients, Cg(r) = 5*^||(y), D{ J(r) = S^(r), and Eg(r) = S% ± (r), where 
we have taken l^ 0 = 1. In effect, these coefficients represent the spherical average over 
the anisotropic particle field. 

Figure 10 shows the probability density function (PDF) for the parallel component of 
the relative velocity for 0 ^ r/rj ^ 2, at different Stokes numbers and a fixed Reynolds 
number of 398. The panel on the left shows the case of no gravity and the panel on the 
right the case of Fr = 0.052. It is immediately apparent that gravity has a profound effect 
on the PDFs, causing a dramatic reduction in the tails of the distribution, particularly 
for particles with larger Stokes numbers. This is consistent with the explanation given in 
§4.1.1, which argued that gravity suppresses path-history effects, thereby reducing the 
relative velocities. 

We next consider the relative velocity variances as a function of the separation distance. 
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Figure 10. PDFs of the radial relative velocities without gravity (a) and with gravity 
(.Fr = 0.052) (b) at Rx = 398 and 0 ^ r/rj ^ 2. The solid line denotes data for fluid (St = 0) 
particles. 



Figure 11. The variance of the particle relative velocities parallel to the separation vector (a) 
and perpendicular to the separation vector (b) for R\ = 398, plotted as a function of r/rj for 
different St. The thin solid lines indicate data with gravity (Fr = 0.052), the dashed lines 
indicate data without gravity, and the thick dotted line indicates fluid (St = 0) particles. The 
Stokes numbers are indicated by the line labels, and the fluid velocity scalings are indicated by 
thick solid lines. The main plots show St = 0 and St is 1, while the insets show St = 0 and 
St = 0.3. 


Figure 11 shows the spherically averaged quantities S^dr) an< ^ ^llM- For St ^ 1 
(corresponding to Sv > 19.3), we observe that gravity strongly decreases the relative 
velocities, by orders of magnitude in some cases. The explanation is that gravity reduces 
the effect of the path-history interactions, causing a reduction in the relative velocities 
(cf. §4.1.1). 

Data at St = 0.3 is shown in the insets in figure 11. For St = 0.3 and very small 
separations (r/ry < 1), the relative velocities show evidence of path-history effects (see 
also Part I), causing the relative velocities to exceed those of St = 0 particles in both the 
longitudinal and transverse directions. However, at larger separations (1 < r/rj < 10), we 
observe a crossover to the expected increase (decrease) in the longitudinal (transverse) 
relative velocities as a result of preferential sampling. We also see that gravity causes the 
relative velocities to approach those of the underlying fluid, since it decreases preferential 
sampling effects. These trends are consistent with the discussion in §4.1.1. [Note that 
since particles with St < 0.3 have only weak gravitational forces (Sv < 5.79), the effect 





19 


Reynolds-number effects on inertial particle dynamics, Part II. 



Figure 12. The mean inward particle relative velocity for R\ = 398, plotted as a function of 
r/ii for different St. The thin solid lines indicate data with gravity (Fr = 0.052), the dashed 
lines indicate data without gravity, and the thick dotted line indicates fluid (St = 0) particles. 
The Stokes numbers are indicated by the line labels, and the fluid velocity scalings are shown 
by thick solid lines. The main plot shows St = 0 and St ^ 1, while the inset shows St = 0 and 
St = 0.3. 

of gravity on the relative velocities is less apparent, and data from these particle classes 
are therefore not shown.] 

The decreased influence of path-history interactions is also apparent in the scaling of 
the relative velocity variances in figure 11. In the dissipation range, the fluid relative 
velocity variances are proportional to r 2 . In the absence of gravity, the relative velocities 
for St ^ 1 do not exhibit r 2 -scaling over any portion of the dissipation range, as was 
noted in Part I. However, with gravity, we see a clearer r 2 -scaling extending much deeper 
into the dissipation range for 1 < r/p < 10. Compared to the case without gravity, one 
has to go to smaller separations to observe deviations from r 2 scaling, as gravity reduces 
the path-history effect responsible for the formation of caustics. 

Of particular interest to the collision kernel is the mean inward velocity (cf. §4.3). 
In general, the mean inward velocity is more difficult to analyze theoretically than the 
relative velocity variance, but qualitatively we expect both statistics to follow the same 
trends. In figure 12, we show S^n(r) at R\ = 398, both with and without gravity. We 
see that the mean inward velocities, like the variances, decrease with the addition of 
gravity at large (small) St due to the reduced influence of the path-history (preferential- 
sampling) mechanisms. Interestingly, with gravity, the mean inward relative velocities (in 
contrast to the relative velocity variances) follow the local fluid scaling almost perfectly 
throughout the entire dissipation range. This result was also noted in Bee et al. (2014). 
One plausible explanation is that since the mean inward velocity is a lower-order statistic 
than the velocity variance, it is less affected by path-history interactions and is more 
affected by the local fluid turbulence, causing the scaling exponents to increase (see Part 

I)- 

To further verify our arguments in §4.1.1, it is helpful to decouple the effects of grav¬ 
ity and inertia by varying each independently. We do so in figure 13, where we show 
S^ii (r/rj = 0.25) and S^fr/r) = 0.25) for 0 < St ^ 56.2, 0 < Sv ^ 100, and R\ = 227. 
While the results at high Sv are likely artificially affected by the periodic boundary con¬ 
ditions (see Appendix A), we can still use these data to discuss the qualitative trends 
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Figure 13. Filled contours of S^Ar/rf)/u v (a) and SLAr/rj)/^ (b) evaluated at r/r/ = 0.25 
and R\ = 227 for different values of St and Sv. The contours are logarithmically scaled, and 
the colorbar labels indicate the exponents of the decade. The diagonal lines denote three differ¬ 
ent values of Fr, corresponding to conditions representative of stratiform clouds (Fr = 0.01), 
cumulus clouds (Fr = 0.05), and cumulonimbus clouds (Fr = 0.3). 


in the relative velocities at different values of St and Sv. We see that S and have 
similar qualitative trends. For Sv < 10, both quantities increase with increasing St, 
either due to preferential-sampling effects (at low St) or path-history effects (at high 
St). The relative velocities also decrease with increasing Sv, since gravity causes both 
effects to be less significant, as discussed in §4.1.1. Finally, we observe that the relative 
velocities are the smallest for the largest values of Sv and St ~ 1. At smaller (larger) 
St, preferential-concentration (path-history) effects are more significant, leading to an 
increase in the relative velocities. An implication of these results for clouds is that for a 
given value of St, droplets in stratiform clouds (Fr = 0.01) will generally have smaller 
relative velocities than droplets in cumulonimbus clouds (Fr = 0.3). 

We have thus far examined and explained relative velocity statistics for fixed Reynolds 
numbers. We now consider how these statistics are affected by changes in R\. In fig¬ 
ure 14, we show both S p _^(r / 1 /) / and S^(r/rj)/u^. For 0 ^ St < 3, the longitudinal 
relative velocities have a very weak dependence on R\, both with and without gravity. In 
Part I, we noted that S^(r/rf) /increases weakly with increasing R\ in the absence of 
gravity for 0.3 < St < 1, and attributed this trend to the increased scale separation with 
increasing R\, which causes a few inertial particles to be affected by their memory of 
increasingly energetic turbulence. (We also suggested that increased intermittency effects 
at higher Reynolds numbers could contribute to the observed trends.) While we expected 
the increased scale separation of the turbulence to also increase the relative velocities for 
higher-Sf particles, we instead found that the relative velocities decreased with increasing 
R\ for 1 < St <3. We argued that this trend was caused by a corresponding decrease in 
the rotation timescales T^ k /t v , which in turn reduced the influence of path-history ef¬ 
fects and decreased the relative velocities. However, with gravity, we find that is 

generally independent of R\ (see §3.1), and thus we expect that this allows the increased 
scale separation of the turbulence to cause S^(r/ij)/u^ to uniformly increase with in¬ 
creasing R\ for St > 0.3. Our results in figure 14 confirm this expectation. We also see 
that S p _^(r/if)/u v is generally unaffected by changes in R\, since it is less influenced by 
the relatively infrequent occurrences of the effects discussed above. 

We are also interested in the effect of the Reynolds number on the transverse relative 
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Figure 14. The normalized mean inward radial relative velocities (a) and relative velocity 
variances (b) plotted as a function of St for different values of R\. The open symbols denote the 
case with gravity (Fr = 0.052), and the filled symbols denote the case without gravity. Data 
are shown for particles with separations r = 0.25r; and r = 9.75 p. 



Figure 15. (a) The transverse relative velocity variances, plotted as a function of St for different 
values of R\. Data are shown for particles with separations r = 0.25p and r = 9.75rj. (b) The ratio 
between the transverse and longitudinal relative velocity variances S% ± (r/rf)/S %y (r/rj), evaluated 
at r/rj = 0.25. In both plots, the open symbols denote the case with gravity (Fr = 0.052), and 
the filled symbols denote the case without gravity. 

velocities, since these statistics will also affect the degree of particle clustering (see §4.2). 
Figure 15(a) shows the transverse relative velocity variances S^^fr/rj) for different values 
of St and R\. We see that the trends with R\ in the transverse direction are identical to 
those in the longitudinal direction. 

In Part I, we found that without gravity, the longitudinal and transverse relative 
velocities became equivalent at small separations for St > 0.3. The physical explanation 
is that the path-history contribution to their relative velocities decreases the coherence 
of the pair motion, and in the ballistic limit where the pairs move independently of each 
other, the longitudinal and transverse components are equal. However, with gravity, path- 
history effects are weaker, and thus the longitudinal and transverse velocities may not 
be the same in this regime. 

We compare these two quantities in figure 15(b) by plotting S^^fr/rj)/S^(r/rj) at 
r/r] = 0.25. We see that both with and without gravity, the ratio S^ir/p)/S^(r/rj) 
approaches 2 at low Stokes numbers, the value for fluid particles (e.g., see Pope 2000). 
For St > 0.3 without gravity, S^^fr/rj) and S^(r/r]) are equivalent, as expected. At 
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Figure 16. The scaling exponents of the longitudinal mean inward relative velocities (a) and 
relative velocity variances (b) plotted as a function of St for different values of R\. The open 
symbols denote the case with gravity (Fr = 0.052), and the filled symbols denote the case 
without gravity. The scaling exponents are computed using linear least-squares regression over 
the range 0.75 ^ r/rj ^ 2.75. 


high values of St with gravity, however, the longitudinal and transverse components are 
not equivalent, since path-history effects are weakened by gravity. As R\ increases, the 
particle relative velocities are affected by increasingly energetic turbulence along their 
path histories. As a result, the relative velocities are larger and the particles move more 
ballistically, causing the ratio S^ir/r])/S^(r/rj) to decrease with increasing R\. 

We now examine the effects of Reynolds number on the scaling of the relative velocities. 
As in Part I, we compute the scaling exponents of the mean inward relative velocity (Cy”) 
and the relative velocity variance (£|) by performing linear least-squares power-law fits 
over separations 0.75 ^ r/r/ ^ 2.75. While the results in figure 11 suggest that the power- 
law exponent may vary considerably over this range, we are forced to use this relatively 
large range to gather sufficient statistics at the smaller separations. Thus, while these 
scaling exponents allow us to assess the trends with R \, they only provide a qualitative 
understanding, as the scaling varies considerably throughout the dissipation range, and 
this variation is not completely captured by this analysis. 

The scaling exponents are shown in figure 16. The trends at a given value of R\ are in 
agreement with our observations from figure 11 and figure 12 , which are explained above. 
We also see that with gravity, the scaling exponents tend to decrease with increasing 
R \, since the particles are more influenced by their memory of increasingly energetic 
turbulence, and therefore tend to move ballistically, as noted above. 

Angular distributions of S ^_,, (r)/S'^|| (r) are presented in figure 17 for a range of Stokes 
numbers, R\ = 398, and r < 77 . The anisotropy of the distribution is reflected in the 
color variation over the surface of the sphere. We see immediately that the asymmetry in 
the relative velocities follows opposite trends at small and large St. That is, at small St, 
the relative velocities are largest for particles separated in the vertical direction, whereas 
for St > 1, the opposite is true. The explanation is as follows. At small separations 
and St < 0.3, particle pairs that are separated vertically will have longitudinal relative 
velocities that are proportional to S 33 (x p (t),t), the longitudinal velocity gradient in 
the vertical direction as sampled by inertial particles (see §3.1). Particle pairs that are 
separated along the horizontal direction, however, will have longitudinal relative velocities 
that are proportional to <Sn(:c p (f), t), the longitudinal velocity gradient in the horizontal 
direction. In §3.1, we observed that particles tend to preferentially sample flow where 
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Figure 17. The directionally dependent mean inward relative velocity ( r ), normalized by the 

spherically averaged mean inward relative velocity Sf_„ (r), shown on a unit sphere for R\ = 398 
and r < y with gravity (Fr = 0.052). The different columns correspond to different values of 
St. The top row shows the projection where gravity is directed into the page, and the bottom 
row shows the projection where gravity is directed downward. 



Figure 18. The second (a) and fourth (b) spherical harmonic coefficients of SO^r), normalized 
by the zeroth spherical harmonic coefficient, plotted as a function of r/y for different St and 
R\ with gravity (Fr = 0.052). The different values of St considered (0.3, 0.5, 1, 3) are shown in 
black, red, blue and green, respectively, and the Stokes numbers are indicated by the line labels. 


the vertical (horizontal) velocity gradients are larger (smaller). As a result, the relative 
velocities are expected to be larger (smaller) for particles which are separated in the 
vertical (horizontal) directions. At large St, however, the relative velocities are smallest 
for particles that are separated vertically. The physical explanation is that the correlation 
timescales of the fluid velocity gradients along this direction are the smallest (see §3.1), 
and thus the particles have a weaker response to vertical fluid velocity fluctuations and 
correspondingly lower particle relative velocities along this direction. 

Next, we consider the dependence of the anisotropy on the separation distance by 
plotting the spherical harmonic coefficients of order 2 and 4 as a function of r in figure 18. 
(Coefficients above order 4 are too small to measure accurately and thus are not shown.) 
We see that the anisotropy in the relative velocities generally decreases with increasing 
separation. The physical explanation is that with increasing separation distance, the 
particle motions are increasingly dependent on larger eddies, and the relative velocities 
induced by these isotropic eddies are increasingly energetic relative to the anisotropic 
velocities induced by gravity. 

For small (large) values of St, c° and c° are positive (negative), indicating that the 
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particle relative velocities are strongest for particles that are separated vertically (hori¬ 
zontally). These observations are in agreement with the trends shown in figure 17. We 
also see that c° and c 3 tend to become more (less) isotropic at high (low) St as R\ 
increases. While the physical explanation for the trend at low St is unclear, at high St, 
we expect that this increase in isotropy is linked to the increase in the relative velocities 
with increasing R\. That is, as the overall relative velocities increase, the anisotropic 
velocities induced by gravity will be comparatively weaker, and thus the particle relative 
velocities will be more isotropic. 


4.2. Particle clustering 

In this section, we consider how gravity affects the inertial particle clustering process. 
We first provide a theoretical explanation for the clustering (§4.2.1) and then compare 
the predictions to DNS results (§4.2.2). 


4.2.1. Theoretical framework for particle clustering 

We use the angular distribution function g(r) (hereafter ‘ADF,’ see Gualtieri et al. 
2009) to quantify the degree and orientation of particle clustering. We define g(r) as 


N(r,9,tj>)/V(r,9,<j>) 

9(0 =- m - 


(4.9) 


In this equation, N(r, 9, <f>) denotes the number of particle pairs in a truncated spherical 
cone with nominal radius r, polar angle 6 , and azimuthal angle f>. The volume of the 
truncated spherical cone V (r, 9, <fi) is given by 

V (r, 9 , 4>) = sm(9)A9A(j) [(r + Ar) 3 — (r — Ar) 3 ] /3, 


where Ar is the radial width, A 9 is the extent of the polar angle, and A (f> is the extent 
of the azimuthal angle. 

As discussed in §4.1.1, the equation governing the relative motion of like particles, (4.1), 
is the same with and without gravity. Hence the theory of Zaichik & Alipchenkov (2009), 
developed without consideration of gravitational settling, is not explicitly changed by the 
presence of gravity; instead gravity modifies the equation through implicit changes to its 
coefficients. From Zaichik & Alipchenkov (2009) (see also Bragg & Collins 2014a), the 
equation describing g(r) at steady-state is 


0 = -StT v (s% + A) • V r g - Str v gV r ■ S%, (4.10) 

where A is the dispersion tensor describing the influence of the fluid velocity difference 
field on the dispersion of the particles (see Bragg & Collins 2014a). 

Since S% is anisotropic with gravity (cf. §4.1.2), so too is the ADF g(r), as is evident 
in (4.10). We similarly adopt an expansion of the ADF in terms of spherical harmonic 
functions 

OO 

ff (r)=£<&(r)l$(0), (4.11) 

£=0 

where C 2(: are the spherical harmonic coefficients for the ADF. We only attempt to model 
the zeroth-order term, which we denote as 


g(r)=CS(r). (4.12) 

where Tq 0 = 1. Note that g(r) is just the spherical average of the ADF g(r), or the 
equivalent of the radial distribution function (RDF) defined in Part I. While this approach 
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captures the leading-order contribution of gravity to the ADF, it does not predict the 
gravity-induced anisotropy. (Note that a similar approach was adopted by Alipclienkov 
& Beketov (2013) to model particle clustering in a homogeneous turbulent shear flow.) 
We will assess the accuracy of this approach in §4.2.2. 

We therefore consider the isotropic form of (4.10), which is given as 


0 = -St 


(S*|l+A||) Vtg-St (v^h+2 r 


qP _ qP 
° 2 \\ ° 2 ± 


(4.13) 


=£> 


where Z denotes a variable Z normalized by Kolmogorov scales. Ay denotes the projection 
of the dispersion tensor A along a direction parallel to the particle separation vector. We 
note that the first term on the right-hand-side is associated with outward diffusion (which 
acts to reduce the RDF), while the second is associated with an inward drift (which acts 
to increase the RDF). The corresponding diffusion and drift coefficients are denoted as 
53 and 5, respectively (see Bragg & Collins 2014a). 

As explained in Bragg & Collins (2014a), the overall degree of clustering is determined 
by the relative strength of the drift and diffusion coefficients, 0/53. Increases (decreases) 
in this ratio are associated with increased (decreased) clustering. To understand how 
gravity affects the clustering, we therefore need to understand its effect upon the ratio 
0/53. ' 

Bragg & Collins (2014a) demonstrated that in the limit St <C 1, 53 and 0 simplify to 

53 « - f 2 B nl , (4.14) 

and 

0 « -Str ((S 2 ) p - (I? 2 ) p ) /3, (4.15) 

respectively, where B n i is the non-local coefficient defined in Chun et al. (2005) that is 
independent of the particle Stokes number. As is evident in (4.14), the diffusion term 
is independent of the Stokes number-hence 53 is unaffected by gravity in this limit. In 
contrast, 0 is proportional to both St and (S 2 ) p — (1Z 2 ) P . Gravity reduces the correlation 
between the particles and the strain and rotation fields, and thereby reduces 0. The 
net effect of gravity in this limit is therefore to reduce the ratio 5/53 and the degree of 
clustering. We also note that at low St, the trends in the strain and rotation rates with 
R\ are very weak. Consistent with Part I, these trends imply the clustering to be nearly 
independent of R\ in the limit St <C 1, with or without gravity. 

For inertial particles with intermediate Stokes numbers in the range 0.2 < St < 0.7 
(note that this range is approximate and is likely sensitive to changes in Reynolds number, 
as mentioned in Part I), we no longer have simple relationships for the drift and diffusion 
terms. However, we can still make qualitative predictions for the trends in the clustering 
based on our physical understanding and the explanations in Bragg & Collins (2014 a). 
Over this range of St, Bragg & Collins (2014a) argued that both preferential-sampling 
effects and path-history effects act to increase clustering, with the path-history effects 
becoming increasingly dominant as St increases. Since gravity decreases both preferential- 
sampling and path-history effects, we therefore expect it to reduce clustering over the 
approximate range of 0.2 < St < 0.7. 

We next use (4.13) to understand clustering at larger values of St. In Part I, we 
argued that without gravity, Sf « S%± for St > 0.3, and we were thus able to neglect 

the term 2r~ 1 (S^ — <Sfj_) in (4.13) at high St. However, with gravity, the longitudinal and 
transverse relative velocity variances are generally not equal (as discussed in §4.1.2), and 



26 


P. J. Ireland et al. 


hence this term must be retained. We are still able to neglect Ay at small separations and 
high St, since this term is inversely proportional to St and decreases as the timescales 
of the fluid velocity seen by the particles decrease (see Bragg & Collins 2014a). 

The simplified form of (4.13) at high St is therefore 


0 = -StS% n Vf9 -St + 2 f~ 


qP _ qP 
° 2 \\ ° 2 ± 


9 • 


—V- 

D 


Taking the ratio between the drift and diffusion coefficients gives us 

5 vA+2f- 1 


qP _ qP 

J 2|l ° 2 ± 


D 


CP 

' 5 2 || 


(Cf+2-2 S p 2± /S p n ), 


(4.16) 


(4.17) 


where (| is the scaling exponent of the longitudinal relative velocity variance (see §4.1.2). 

From §4.1.2, we see that with gravity, increases very strongly, while ‘IS^^/S^ in¬ 
creases weakly. We expect that £| — ZS^/S^ will therefore increase with gravity, leading 
to an increase in the clustering at high St. We attribute these trends to the fact that 
gravity reduces the influence of path-history effects on particle relative motion in the 
dissipation range at high St. Path-history effects primarily act to diminish clustering 
for St > 1 (see Bragg & Collins 2014a), and thus gravity reduces their influence and 
therefore acts to increase g{f). We emphasize that this does not necessarily imply that 
gravity increases the inward drift and decreases the outward diffusion. In fact, our data 
indicate that gravity decreases both 5 and 2). The reduction in 2) with gravity, however, 
is stronger than the corresponding reduction in 5, causing the ratio D/2), and therefore 
the clustering strength, to increase. We also see from §4.1.2 that increasing R\, with 
gravity, generally leads to a decrease in both and S/V/S/,/, due to the increased in¬ 
fluence of path-history effects on particle relative motion at higher Reynolds numbers. 
While both quantities seem to be decreased by about the same amount with increasing 
R \, the latter quantity has a greater effect on the clustering, since it is multiplied by a 
factor of two in (4.17). This explains why 5/2) decreases as R\ increases, leading to a 
decrease in the clustering. 

In summary, gravity reduces clustering at low and intermediate values of St, but in¬ 
creases clustering at high values of St. Gravity reduces both preferential-sampling and 
path-history effects. Preferential-sampling effects increase clustering at low values of St, 
and thus gravity reduces clustering by weakening the preferential sampling drift mecha¬ 
nism. At intermediate values of St, both preferential sampling and path-history effects 
act to increase clustering, and thus gravity, by weakening both effects, reduces cluster¬ 
ing. At higher values of St, path-history effects act to diminish clustering, and gravity 
increases clustering by decreasing path-history effects. 

This explanation is consistent with arguments put forth by Franklin et al. (2007); 
Ayala et al. (2008); Woittiez et al. (2009); Onishi et al. (2009); Rosa et al. (2013) for 
low-Stokes-number particles. However, at higher Stokes numbers, several of these authors 
have argued that gravity facilitates interactions with large-scale turbulent eddies, which 
in turn leads to increased particle clustering (Franklin et al. 2007; Woittiez et al. 2009; 
Rosa et al. 2013). In contrast, we argue that gravity reduces the temporal correlation 
radius over which the particles are affected by their path-history interactions with the 
turbulence, and therefore causes the particle relative motion in the dissipation range to 
be less affected by their interaction with larger-scale turbulence. 

Another recent paper by Park & Lee (2014) suggests the increase in particle clustering 
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Figure 19. Instantaneous particle locations in 10007/ x 1000?/ x 1077 slices of the domain for 
different values of St and R\ = 398. The top row is for the case without gravity, and the bottom 
row is for the case with gravity (Fr = 0.052), with the gravitational vector pointing downward. 


with gravity at high Stokes numbers is linked to the skewness of the vertical velocity gra¬ 
dients of the underlying fluid. This proposed explanation, however, seems to be flawed, 
and we suggest that their finding reflects a certain correlation in the system, but not 
a causal relation. This conclusion is supported by the recent work of Gustavsson et al. 
(2014), who observed qualitatively similarly strong clustering of high-Stokes-number par¬ 
ticles in the presence of gravity in their random Gaussian flow field, whose fluid velocity, 
by definition, has no skewness. 

4.2.2. Particle clustering results 

We begin by first examining instantaneous snapshots of particle positions in the sim¬ 
ulations in figure 19. From these visualizations, it is evident that gravity alters both the 
degree and orientation of the clusters, and generally causes them to be aligned with the 
gravity vector. Recall that the clustering observed at St > 1 is not related to the Maxey 
centrifuge mechanism (Maxey 1987), but rather is due to the history effect (Bragg & 
Collins 2014a). This explains why liigher-Stokes-number particles show signs of cluster¬ 
ing even though they nearly uniformly sample the strain and rotation fields (cf. figure 2). 
We conclude the vertical streaks in figure 19 are a manifestation of the non-local history 
mechanism acting in the presence of gravity. 

To quantify the clustering, we first consider the spherically averaged RDFs g(r) to 
analyze the overall effect of gravity on the degree of clustering. Figure 20 shows plots of 
the RDFs for St = 0.3 and St = 3 both with gravity ( Fr = 0.052) and without gravity 
at R\ = 398. The results show gravity causes the RDFs to decrease (increase) at low 
(high) values of St, in agreement with our arguments in §4.2.1 and the earlier findings 
(Franklin et al. 2007; Ayala et al. 2008; Woittiez et al. 2009; Onishi et al. 2009; Rosa 
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Figure 20. DNS (symbols) and theoretical (lines) data for the RDFs g plotted as a function of 
r/r) at R\ = 398 both with gravity (open symbols and solid lines) and without gravity (filled 
symbols and dashed lines), (a) shows RDFs for St = 0.3, and (b) for St = 3. The theoretical 
predictions are calculated from the equations in Zaichik & Alipchenkov (2009) with the non-local 
correction from Bragg & Collins (2014a), using our DNS data to specify T$ s , Sl/^, and S^x- 

et al. 2013; Bee et al. 2014; Gustavsson et al. 2014). For sufficiently large separations, 
the RDFs with and without gravity approach each other. 

To test the theory in §4.2.1, we compute g(r) from (4.13) using the relative velocity 
statistics 5',^ and S% ± obtained from DNS. In addition, we use the directionally averaged 
strain timescales Tg S from the DNS (see §3.1) and the non-local closure proposed in 
Bragg & Collins (2014a) to compute the dispersion tensor. This allows us to test the 
formulation of (4.13) and the theoretical arguments in §4.2.1. Figure 20 shows that the 
theory captures the quantitative results in the DNS well, indicating that (4.13) is an 
accurate model even for an anisotropic particle phase, and thus verifying the physical 
explanations presented in §4.2.1. 

Next, we consider the dependence of the RDFs on R\ in figure 21. As was the case 
without gravity (see Part I), the RDFs are largely independent of R\ for St < 1 with 
gravity (figure 21(a)). This is in agreement with our arguments in §4.2.1 and implies that 
small-S't clustering is a small-scale phenomenon (both with and without gravity) that is 
generally unaffected by the intermittency of the turbulence. At larger St (figure 21(b)) 
the RDFs increase monotonically with increasing R\, since the ratio between the drift 
and diffusion increases, as discussed in §4.2.1. 

We can also quantify the degree of small-scale clustering by performing a power-law 
fit of the RDFs (Reade & Collins 2000) at small separations, 

g{r) « c 0 0) , (4.18) 

as discussed in Part I. The power-law fits are performed over the range 0.75 ^ r/r] ^ 2.75, 
and the calculated values of Co and c\ are plotted in figure 22. We observe that both the 
power-law coefficient Co and the exponent c\ decrease when gravity is introduced for 
St < 1.5 and increase when gravity is introduced for St > 1.5, consistent with our 
explanations in §4.2.1. The theoretical model for g(r/r]) from (4.13) (with the relative 
velocities and strain timescales specified from the DNS) is in good agreement with the 
DNS data, both with and without gravity. We also note that our DNS results for the 
exponent c\ agree well with those of Rosa et al. (2013) and Bee et al. (2014) (not shown). 

To further investigate the two-parameter space of particle inertia and gravity, we con- 
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Figure 21. The RDFs g(r/rf) with gravity (Fr — 0.052) for different R\ for St ^ 1 (a) and for 
St — 3 (b). The Stokes numbers in (a) are indicated by the line labels. 




Figure 22. The prefactor Co (a) and the exponent ci (b) of the power-law fits of the RDFs 
using (4.18). Open symbols denote data with gravity (Fr = 0.052), and filled symbols denote 
data without gravity. The predictions from Zaichik & Alipchenkov (2009) for R\ = 398 (where 
DNS data are used to specify the relative velocities and the strain timescales here) are shown 
with solid lines (gravity) and dashed lines (no gravity). 


sider the RDFs at R\ = 227 for 0 < St ^ 56.2 and 0 < Sv f 100 in figure 23. While 
the particle clustering behavior here is generally complex and varies strongly with St 
and Sv, we are able to provide physical explanations for several of the observed trends. 
For St ^ 1, the RDFs generally decrease with increasing St, since the particles become 
unresponsive to almost all of the underlying turbulence, as expected. Also, for St < 1, 
the RDFs tend to decrease with increasing Sv, since the preferential-sampling and path- 
history mechanisms may act to increase clustering here, and they are both reduced by 
gravity. We also see that for 1 < St < 10, the RDFs generally increase with increasing 
Sv, as expected, since gravity increases the ratio between the drift and the diffusion, as 
discussed in §4.2.1. 

We now investigate the anisotropy of the particle field by introducing higher-order 
spherical harmonic functions, 


g( r ) 

g( r ) 


OO 
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c 2e( r ) 

C 0 °(r) 




in figure 24 for various values of St with gravity (Fr = 0.052) at R\ = 398. 


(4.19) 
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Figure 23. Filled contours of the RDF evaluated at r/rj = 0.25 (a) and r/rj = 9.75 (b) for 
different values of St and Sv, for R\ = 227. The diagonal lines denote three different values of 
Fr, corresponding to conditions representative of stratiform clouds (Fr = 0.01), cumulus clouds 
(Fr = 0.05), and cumulonimbus clouds (Fr = 0.3). 


St=0.1 St=0.3 St=0.5 St=0.7 St=l St=3 



Figure 24. The ADF g ( r ) (normalized by the RDF g(r )) shown on a unit sphere for Rx = 398 
and separations r < g with gravity (Fr = 0.052). The different columns correspond to different 
values of St. The top row shows the projection where gravity is directed into the page, and the 
bottom row shows the projection where gravity is directed downward. 


Consistent with the qualitative observations from figure 19, we see that particle cluster¬ 
ing is strongest along the vertical direction, in agreement with earlier findings of Dejoan 
& Monchaux (2013), Bee et al. (2014) and Park & Lee (2014). At low St, the anisotropic 
clustering is caused by the fact that particles tend to preferentially sample downward- 
moving flow (see §3.2). Bee et al. (2014) showed that this preferential sampling causes 
particles to form vertical clusters when St is small. 

When St is large, the effects of preferential sampling vanish, and the anisotropy is 
related to the way gravity influences the path-history effect. Since the analysis presented 
in §4.2.1 is restricted to spherically averaged clustering, the results cannot be used to 
predict the trends in anisotropy. However, we hope that the strain and rotation timescale 
results presented in §3.1 will aid in the development of future models to predict this 
anisotropy. 

We also note that the anisotropy increases with increasing St. The physical explanation 
is that as St (and thus Sv) increases, gravitational forces become more significant, causing 
the particle motion to be more anisotropic. However, we expect that at some sufficiently 
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Figure 25. The second (a) and fourth (b) spherical harmonic coefficients of the angular distri¬ 
bution function, normalized by the zeroth spherical harmonic coefficient, plotted as a function 
of r/r), for different St and R\ with gravity (Fr = 0.052). The different values of St considered 
(0.3, 0.5, 1, 3) are shown in black, red, blue and green, respectively, and the Stokes numbers are 
indicated by the line labels. 


large value of St and Sv, the particles are unaffected by the fluid turbulence, and thus 
clustering (and the associated anisotropy) will vanish. 

We explore the degree of anisotropy as a function of the separation distance r and the 
Reynolds number by plotting the spherical harmonic coefficients C® and in figure 25 
(coefficients above order four are too small to be statistically significant, and hence are not 
shown). In agreement with figure 19, we see that the anisotropy increases with increasing 
St. We also observe that as r/r/ increases, the anisotropy approaches zero, since both the 
clustering and the clustering anisotropy vanish in the limit r/r] —> oo. 

The degree of anisotropy is also Reynolds-number-dependent and decreases with in¬ 
creasing R\ for larger values of St. Since the relative velocities of high-S't particles gen¬ 
erally also tend to become more isotropic as R\ increases (see §4.1.2), it is likely that 
the reduction in the anisotropy of the relative velocities will cause a similar reduction in 
the anisotropy of the clustering. 

4.3. Particle collision kernels 

The final two-particle statistic we consider is the kinematic collision kernel K for inertial 
particles. Sundaram & Collins (1997) and Wang et al. (2000) showed that for an isotropic 
particle field, K is given by 

K{d) = 4tt <d 2 g(d)S p _ ll (d), (4.20) 

where d is the particle diameter. In Appendix C, we show mathematically that (4.20) 
holds even for an anisotropic particle field. As in Part I, we plot the non-dimensional 
collision kernel K{d) = K(d)/(d 2 u ri ) in figure 26 (see also Vofikuhle et al. 2014). Note 
that while we simulate only point particles, we define d from St for this plot by prescrib¬ 
ing the density ratio to be p p /pf = 1000. Since our statistics are generally not sufficient 
to compute the RDFs and relative velocities at separations on the order of the parti¬ 
cle diameter, we fit both quantities using linear least-squares power-law regression and 
extrapolate the resulting power-law fits to r = d (refer to Rosa et al. 2013). 

Figure 26 shows the dimensionless collision kernel, with and without gravity, as a 
function of the Stokes number. Notice that the collision kernel is reduced by gravity. At 
low St, the RDF and the relative velocity both decrease with gravity, thereby decreasing 
the collision kernel. At higher St, the reduction in the relative velocity with gravity (see 
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Figure 26. The non-dimensional collision kernel K(d) plotted as a function of St for 
pp/pf = 1000 and different values of R\. The open symbols denote data with gravity 
(Fr = 0.052), and the filled symbols denote data without gravity. 



Figure 27. The non-dimensional collision kernel K(d) plotted as a function of St for different 
values of p v /pf and R\. All data are with gravity. Filled red symbols indicate p P /pj = 250, 
open symbols indicate p p /pf = 1000, and filled blue symbols indicate p P /pf = 4000. 


§4.1.2) is sufficient to compensate for the slight increase in the RDF with gravity (see 
§4.2.2), sustaining the net reduction in the collision kernel. Finally, we note that at large 
St with gravity, the collision rates increase with increasing R\ since both the relative 
velocity and RDF increase with increasing Reynolds number. 

We found in Part I that K(d) without gravity is nearly independent of p P /pf for 
St > 1, in agreement with Vofikuhle et al. (2014). The physical explanation is that g(d) 
and S'0||(<i)/tt ?; are either independent of d (for St ^ 10) or have inverse power-law 
scalings with d that precisely cancel (for 1 < St ^ 3). In contrast, with gravity, path- 
history effects are suppressed, and these quantities have different scaling behaviors (see 
§4.1.2 and §4.2.2). We therefore find that with gravity, A (d) decreases with increasing 
Pp/pf at all values of St considered, as shown in figure 27. 
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Figure 28. Filled contours of K(d) for different values of St and Sv at R\ = 227 for 
pp/pf = 1000. The contours are logarithmically scaled, and the colorbar labels indicate the 
exponents of the decade. The diagonal lines denote three different values of Fr, corresponding 
to conditions representative of stratiform clouds (Fr = 0.01), cumulus clouds (Fr = 0.05), and 
cumulonimbus clouds (Fr = 0.3). 

Next, we expand our parameter space to consider the collision kernel for 0 < St ^ 56.2 
and 0 < Sv ^ 100 at R\ = 227 in figure 28. Here, we limit our analysis to cases where 
Pp/pf = 1000. Once again, we caution that the statistics at the highest values of St and 
Sv may be affected by the periodic boundary conditions in the vertical direction (see 
Appendix A), the use of a linear drag model (cf. §2.2), and the extrapolations based on 
power-law fits that become questionable for St > 3. Nevertheless, the contour plot shown 
in figure 28 captures the qualitative trends of the collision kernel. In agreement with the 
discussion above, we see that the collision kernel generally increases with increasing St 
and decreases with increasing gravity. A practical implication is that at a given value 
of St, cumulonimbus clouds (Fr « 0.3) will have more frequent droplet collisions and 
growth than stratiform clouds (Fr ss 0.01). 

We now consider the anisotropic collision kernel, K(r, 6, <f >), defined as 

K(d, 6, </>) = ±nr 2 g(d, 9, 0)^„ (d, 9, cf). (4.21) 

K(d, 9, cf>) provides a measure of the rate at which particles with diameter d collide along 
an orientation defined by the radial angle 9 and azimuthal angle cf. We plot the ratio 
between K(d,9,<j>) and its spherical average in figure 29 at R\ = 398. (We are unable 
to show the anisotropic collision kernel at r = d due to inadequate statistics at these 
small separations. We instead show data at 0 ^ r < rj, the smallest separation range over 
which adequate statistics are available.) 

At low St, the collision kernel is approximately isotropic, since gravitational effects 
are weak. Interestingly, the collision kernel also tends toward isotropy at large St, due 
to the opposing trends in anisotropy of the relative velocities (see figure 17) and ADFs 
(see figure 24). For St between 0.1 and 1, particles are more likely to collide along the 
vertical direction, since both the ADFs and the relative velocities are strongest in the 
vertical direction for these values of St. 

Before closing this section, we emphasize two practical implications of these collision 
results for the cloud physics community. The first is that since the DNS indicates that the 
collision rates of particles with low and moderate St are independent of R\, it is likely 
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Figure 29. The anisotropic collision kernel I\(r,9,(f)) (normalized by the spherical average 
K(r)) shown on a unit sphere for R\ = 398 and separations r < r/ with gravity (Fr = 0.052). 
The different columns correspond to different values of St. The top row shows the projection 
where gravity is directed into the page, and the bottom row shows the projection where gravity 
is directed downward. 


that the collision kernels computed here will be useful for predicting droplet collisions 
in high-Reynolds-number atmospheric clouds. The second is that gravity significantly 
reduces the collision kernels for St > 0.1, which implies that simulations without gravity 
can significantly over-predict the collision rates of droplets in atmospheric clouds, and 
highlights the need to include gravity in the analysis of collisional droplet growth in 
turbulent clouds. 


5. Conclusions 

In this study, we explored the influence of gravity on inertial particle statistics in 
isotropic turbulence, with the goal of understanding the turbulence mechanisms con¬ 
tributing to droplet collisions. Our simulations were performed over the largest Reynolds- 
number range (90 < R\ < 597), domain lengths {£/£ < 40), and range of particle classes 
(0 < St ^ 56.2, 0 < Sv < 100) to date. We showed that such large domain sizes are 
essential to obtain accurate statistics of settling inertial particles, suggesting that earlier 
published DNS studies may have errors due to the periodic boundary conditions in the 
vertical direction (see Appendix A). 

Our results indicate that preferential sampling affects the dynamics of particles with 
St -C 1, both with and without gravity. Gravity decreases the degree of preferential 
sampling by reducing the time the particles have to interact with the underlying turbulent 
fields. In particular, gravity reduces the Lagrangian timescales for strain and rotation 
along particle trajectories. As gravitational forces increase, the particles fall more rapidly 
through the flow, leading to smaller velocities and larger accelerations than in the case 
without gravity. We developed models for the Lagrangian strain and rotation timescales 
along particle trajectories and the particle acceleration variances that are valid in the 
limit Sv u'/urj. The model predictions are in reasonable agreement with the DNS in 
this limit. We also find that the mean settling velocity is independent of R\ at low St, 
in agreement with Bee et al. (2014). 

We then applied this understanding of the gravity effect on single-particle statistics to 
two-particle statistics relevant for predicting the collision kernel. At high St, we observed 
that gravity reduces the particle relative velocities from their values without gravity by 
reducing the path-history effect. At low St, gravity acts primarily to reduce the degree 
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of preferential sampling, causing the relative velocities at small separations to be closer 
to those of fluid particles. 

Next, we relate the trends in the relative velocities to those in the particle clustering by 
considering the effect of gravity on a derivative model of the theoretical work of Zaichik 
& Alipchenkov (2009). With gravitational effects included, the model is able to predict 
the spherically averaged RDFs very accurately when DNS data is used to prescribe the 
relative velocities and the Lagrangian strain timescales. By analyzing the model at low 
St, we see that the primary effect of gravity is to decrease the inward drift, leading to 
a decrease in the RDFs. At higher St, in contrast, gravity modulates both the inward 
drift and outward diffusion by diminishing the path-history effect in such a way as to 
increase the drift-to-diffusion ratio, causing an increase in clustering with gravity. We 
also find that the degree of clustering is only weakly dependent on R\ at low St, but 
becomes increasingly sensitive to R\ at higher St. The model captures both of these 
trends through its internal description of the path-history effects. 

We also quantify the degree of anisotropy in these two-particle statistics using a spher¬ 
ical harmonic decomposition. Our results indicate that the particle angular distribution 
functions and radial relative velocities can have anisotropies on the order of 25%, with 
the degree of anisotropy generally peaking at small r/rj. At larger separations, the rela¬ 
tive velocities induced by turbulence become comparatively larger, limiting the effect of 
gravity on particle dynamics. 

We used these data for the RDFs and relative velocities to compute the particle collision 
kernel. As in Part I, we find that the collision kernel is generally independent of R\ at low 
St, while it increases with increasing R\ at high St. We analyze the collision kernel using 
spherical harmonic decompositions, and find that the collision kernel is generally more 
isotropic than the ADFs and the mean inward relative velocities, since the anisotropies 
in the ADFs and the relative velocities have opposing trends at large St. 

We conclude by highlighting some practical implications of this work for the cloud 
physics and turbulence communities and suggesting promising research directions. As 
in Part I, the fact that the collision rates are generally independent of R\ for droplet 
sizes representative of those in warm, cumulus clouds suggests that the collision rates 
predicted here (for R\ < 600) may be representative of those found in atmospheric 
clouds (R\ ~ 10,000). Of course, more simulations at even higher Reynolds numbers 
would help to verify this conjecture. Also, as noted in §4.3, we observed that collision 
rates with gravity are considerably lower than without gravity, suggesting that earlier 
work that neglected gravity may have over-predicted the collision kernels. 

While our study has focused on like-sized particles, our results suggest that gravity will 
have a significant effect on the collision kernels of unlike (different-sized) particles. For 
example, we observe that gravity tends to enhance the settling speeds of particles with low 
and intermediate values of St. The coupling between turbulence and gravity may thereby 
increase the relative velocities of unlike particles, leading to more frequent collisions. In 
addition, we find that gravity causes inertial particles in turbulence to experience large 
accelerations. Chun et al. (2005) found that large accelerations can contribute to higher 
relative velocities between unlike particles, leading to higher collision frequencies. 

Finally, one especially promising finding of this study is that by extending the model 
of Zaichik & Alipchenkov (2009) and using DNS data to specify the Lagrangian strain 
timescales and relative velocities, we are able to accurately predict the RDFs both with¬ 
out gravity (as was also found in Bragg & Collins 2014a) and with gravity. As noted in 
Bragg & Collins (20146), however, this model provides poor predictions of the relative 
velocities of particles with moderate-to-large inertia. Future work should therefore be 
directed at improving these relative velocity predictions. In addition, while we developed 
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a model for the strain and rotation timescales in the limit of large St with gravity, we 
generally require a model for all values of St. If accurate models can be developed for 
the relative velocities and the strain and rotation timescales over the entire range of St 
in atmospheric clouds, we would be able to use the theory of Zaichik & Alipchenkov 
(2009) to model the collision kernel accurately, both with and without gravity, without 
requiring any inputs from DNS. 
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Appendix A. Periodicity effects 

As noted in Woittiez et al. (2009), particle statistics in DNS may be artificially influ¬ 
enced by the periodicity of the domain when gravitational forces are strong. In particular, 
Woittiez et al. (2009) estimated that periodicity effects become significant when the time 
it takes a particle to fall through the domain (~ Cr~ l g~ 1 ) is less than the large-eddy 
turnover time Tj J = £/u', or equivalently, when 


St. Cv!_ 
Fr ~ I u v 


(Al) 


In this case, a particle can artificially encounter the same large eddy multiple times due to 
the finite domain length. We define St cr it as the value of St at which St/Fr = Cu'/(£u v ). 
This roughly corresponds to the largest Stokes number at which we can expect the results 
to be unaffected by periodicity. 

To study the effects of periodicity, we systematically increased the domain size for 
the cases at the three lowest Reynolds numbers, while keeping the large scales, small 
scales, and the forcing parameters the same. The simulation parameters are summarized 
in table 3. For convenience, we will hereafter refer to the simulations from groups I, II, 
and III by their nominal Reynolds numbers of R\ « 90, 147, and 230, respectively. In all 
cases, 0 < St < 3 and Fr = 0.052. 

The resulting energy spectra are shown in figure 30. In all cases, forcing is applied 
deterministically to wavenumbers with magnitude n = -\/2, and thus the location of the 
peak in the spectra is approximately the same for all domain sizes. The spectra are nearly 
identical to the right of the peak, indicating that the inertial- and dissipation-range flow 
scales remain unchanged. We notice that as the domain size increases, larger-scale (i.e., 
lower-wavenumber) flow features become present. These features, however, account for 
very little of the overall energy, and thus they do not significantly affect the fluid statistics. 

We first consider the effect of the domain size on the strain and rotation rates along 
particle trajectories. ( S 2 ) p and (1Z 2 ) P (not shown) are generally invariant with changes 
in the domain size, since they are small-scale, single-time quantities. The Lagrangian 
strain and rotation timescales in figure 31, however, appear to be over-predicted on the 
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Table 3. Simulation parameters for the periodicity study. 

I II III 


R\ 

88 

90 

90 

90 

140 

145 

147 

226 

230 

N 

128 

256 

512 

1024 

256 

512 

1024 

512 

1024 

C 

2tt 

4-7T 

871 

16?r 

2tt 

4-7T 

871 

2tt 

4-7T 

Stcrit 

1.06 

2.06 

4.26 

8.55 

1.39 

2.69 

5.55 

1.78 

3.39 

£ 

1.46 

1.52 

1.47 

1.47 

1.40 

1.49 

1.44 

1.40 

1.49 

t/p 

v! 

55.8 

57.7 

55.7 

55.6 

106 

111 

107 

202 

213 

0.914 

0.912 

0.912 

0.912 

0.914 

0.916 

0.914 

0.915 

0.914 

u'/u v 

4.77 

4.81 

4.82 

4.82 

6.01 

6.12 

6.15 

7.60 

7.70 



10 _1 10 ° 10 1 10 2 10 3 


K 

Figure 30. Energy spectra for different domain lengths C and nominal Reynolds numbers 
R\ = 90 (a), R\ = 147 (b), and R\ = 230 (c). All values are in arbitrary units. 


smaller domain sizes at large St. The explanation is that these timescales involve integrals 
of the autocorrelations of the strain and rotation components (see Equation 3.1), and 
artificially high correlations can result when a particle is wrapped around a domain 
boundary, increasing the timescales. At the highest values of St, our results suggest that 
the strain and rotation timescales have not yet converged to a grid-independent result, 
even though our scaling argument in (Al) indicates that the Stokes numbers considered 
are generally well-below St CI i t . Since these timescales are very small (~ O.lr^), their 
values can be significantly affected by very weak correlations induced by periodicity, and 
thus extremely large domain sizes may be necessary for these artificial correlations to 
vanish. We also see that our DNS timescales are considerably larger than the theoretical 
predictions from §B.l. While part of this discrepancy at large St is due to periodicity 
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st 


Figure 31. The Lagrangian strain (open symbols) and rotation (filled symbols) timescales for 
different domain lengths C and nominal Reynolds numbers R\ = 90 (a), R\ = 147 (b), and 
R\ = 230 (c). The theoretical predictions for the strain and rotation timescales for Sv u' /it,, 
are shown with solid and dashed lines, respectively. 


effects, it is also possible that the theory in §B.2 requires larger values of Sv than we are 
able to simulate here. 

We next discuss the effect of periodicity on velocity and acceleration statistics. The 
turbulence-induced changes in the mean particle settling speeds (not shown) are inde¬ 
pendent of the domain size to within statistical noise. This is presumably because the 
turbulence-induced settling speed modifications of large-iSf particles (which experience 
the strongest periodicity effects) are negligible, and thus the domain periodicity will not 
lead to significant changes in this statistic. It is possible, however, that the settling speeds 
with a nonlinear drag model (which are reduced at large St, see Good et al. 2014) will 
be affected by the domain size. 

We show acceleration variance statistics in figure 32. The acceleration variances are 
under-predicted at the smallest domain sizes, and converge to the theoretical predictions 
at high St (see §B.2) as the domain size is increased. The vertical acceleration variances 
are more affected by the periodic boundary conditions than the horizontal acceleration 
variances, since the vertical velocities are correlated over longer lengthscales (as explained 
in §3.3). 

We conclude this section by examining the relative velocity statistics and the RDFs. 
Figure 33 and figure 34 show the longitudinal mean inward relative velocities and relative 
velocity variances, respectively. These statistics are only weakly sensitive to the domain 
size for all St and converge to a grid-independent result. The relative velocity variances 
appear to be slightly more sensitive to the domain size than the mean inward relative 
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Figure 32. DNS data for {a\) p /a“f (filled symbols) and (a%) v /a^ (open symbols) for different 
domain lengths C and nominal Reynolds numbers R\ = 90 (a), R\ = 147 (b), and R\ = 230 
(c). The theoretical predictions for the horizontal (Equation Bll) and vertical (Equation B12) 
acceleration variances for Sv u'/u v are shown with dashed and solid lines, respectively. 


velocities, presumably because the former are higher-order statistics and are thus sensitive 
to larger and more intermittent flow features. 

Figure 35 indicates that the RDFs are almost entirely unaffected by the finite domain 
sizes at low St. At St = 3 (figure 36), however, the RDFs are evidently quite sensitive to 
the domain size, and do not converge to a grid-independent result. It is possible that these 
RDFs are influenced by the small fraction of larger-scale features which appear as grid 
size is increased (see figure 30), or that the scaling argument in (Al) is not sufficiently 
stringent for the RDF statistics. In any case, the results at large St should be interpreted 
with caution, since the linear drag model and the point-particle approximation are less 
accurate here (see §2.2). 


Appendix B. Theoretical models in the strong-gravity limit 

In this section, we derive theoretical models which are expected to be valid in the 
limit of strong gravitational forces. We require two conditions for these models to hold. 
First, particles must fall so quickly that the spatial decorrelation of the turbulence along 
their trajectories occurs much faster than the Eulerian timescale of the flow. Second, a 
particle’s displacement over a time s is to leading order equal to r p gs. 

The conditions can be shown to be satisfied when Sv v!/u v . It is important to note 
that u '/~ R.f in isotropic turbulence, which implies that our models will only hold 
for very large droplets (i.e., very large values of Sv) in high-Reynolds-number atmospheric 
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Figure 33. The longitudinal mean inward relative velocities (normalized by u v ) for different 
domain lengths L for nominal Reynolds numbers R\ = 90 (a), R\ = 147 (b), and R\ = 230 (c). 
The different Stokes numbers considered (St = 0.3, 1, 2, 3) are shown in black, red, blue, and 
green, respectively, and the Stokes numbers are indicated by the line labels. 


Table 4. Estimated atmospheric-cloud conditions at which our strong-gravity timescale model 

is expected to hold. 


Cloud type Stratiform Cumulus 
e (m 2 /s 3 ) 10“ 3 10“ 2 

R\ 10,800 15,900 

u'/ur, 52.8 64.0 

d (fim) S> 127 >> 186 


Cumulonimbus 

1CT 1 

23,300 

77.6 

» 273 


clouds. In table 4, we provide an estimate of the droplet diameters d at which our model 
is expected to be valid for three different cloud types: stratiform clouds (e ~ 10 -3 m 2 /s 3 ), 
cumulus clouds (e ~ 10~ 2 m 2 /s 3 ), and cumulonimbus clouds (e ~ 10 _1 m 2 /s 3 ). 

From table 4, we see that the droplet diameters at which these models are expected 
to be valid (on the order of hundreds to thousands of microns) are well above those of 
droplets in the ‘size-gap’ in atmospheric clouds (between about 30 and 80 microns-see 
Part I), where turbulence is expected to play a key role in droplet growth. However, while 
the models are not expected to yield accurate predictions of droplet dynamics in the size- 
gap, they provide an important asymptotic limit in which to understand heavy-particle 
motion. Furthermore, since the Reynolds numbers in our simulations are well below those 
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r/v 

Figure 34. The longitudinal relative velocity variances (normalized by iff) for different domain 
lengths C for nominal Reynolds numbers R\ = 90 (a), R\ = 147 (b), and R\ = 230 (c). The 
different Stokes numbers considered (St = 0.3, 1, 2, 3) are shown in black, red, blue, and green, 
respectively, and the Stokes numbers are indicated by the line labels. 


in atmospheric clouds, we expect the models to agree well with our simulation data at 
the largest values of Sv. 

In §B.l, we develop a model for the Lagrangian strain and rotation timescales of the 
flow seen by inertial particles, and in §B.2, we model the particle accelerations. 


B.l. Lagrangian strain and rotation timescales along particle trajectories 
In the limit Sv u r /u v , we can model the Lagrangian strain timescale Tg. . Sk as 

OO 

T SijS km = 7 - 1 -r / [ {s ij (X,s)S km (0,0))e(X,s)dXds, (Bl) 

where g(X , s) = — t p qs). Here and through the rest of this section, no summation 

is implied by repeated indices in the timescale expressions. 

Note that (Bl) is constructed by assuming that the turbulence evaluated along the 
particle trajectories is uncorrelated with the trajectory itself (i.e., ‘Corrsin’s hypothesis,’ 
see Corrsin 1963). In the limit Sv v!/u v , Corrsin’s hypothesis is expected to become 
exact, since the particle motions are almost independent of the underlying flow field here. 
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r/v 

Figure 35. RDFs for different domain lengths C for nominal Reynolds numbers R\ = 90 (a), 
R\ = 147 (b), and R\ = 230 (c). Data are shown for St ^ 1, with the Stokes numbers indicated 
by the line labels. The different Stokes numbers considered (St = 0.1, 0.3, 0.5, 0.7, 1) are shown 
in black, red, blue, green, and magenta, respectively, and the Stokes numbers are indicated by 
the line labels. 


Taking the integral of (Bl), we obtain 

OO 

^SijSkm = 7 \“ [ 0)iSfc m (0,0)\ ds. (B2) 

To derive an analytical expression for the timescales, we assume that ( Sij(r p Qs , 0)<Sfc m (0, 0)) 
decorrelates exponentially, which gives us 


T s ij s km =ts ij s km ,3/(T P Q). (B3) 

£s r j s km ,3 is the integral lengthscale of SijSkm evaluated along the a^-direction. We express 
(B3) in non-dimensional form, 


rpP 

1 SijSkm 


tSijSkr, 


Sv 


(B4) 


where we have used the top-hat symbol to denote a variable normalized by Kolmogorov 
scales. The rotation timescales ^ are defined analogously, 
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Figure 36. RDFs for different domain lengths C for nominal Reynolds numbers R\ = 90 (a), 
R\ = 147 (b), and R\ = 230 (c). Data are shown for St = 3. 


B.2. Inertial particle acceleration variances 

We now develop a model for the particle acceleration variances in the limit Sv v!/u v . 
Following Bee et al. (2006), we write the formal solution for the particle acceleration as 


a p (t) = \ [ exp 

7~p J —OO 


S — t 


[u(x p (t),t) — u(x p (s ), s)] ds. 


(B6) 


Without loss of generality, we take t = 0 and write the particle acceleration covariance 
tensor (a p (0)a p (0)) as 


a p (0)a p (0)) =—j ( u(x p (0),0)u(x p (0),0) 


o 

5 — 

o 

s 


'S' 

/ / exp 

J — OO J — OO 

Jp. 

exp 

T P- 


dsdSj 


-"s 

o 

o 

s 


'S' 

/ / exp 

J —OO J — OO 

Jp. 

exp 

T P. 


u(x p (s ), s)dsdS j 


o 

o 

: — 

r si 


IS] 

/ / exp 

/ — OO J —OO 

1 1 ? 
_1 

exp 

1 ^ 
_i 


u(x p (s), s)u(x p (S), S)dsdS^ . 


(B7) 
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By assuming that preferential concentration effects are weak, we can approximate the 
first term on the right-hand-side of (B7) as 

r- 2 /tx(^(0),0)u(^(0),0)) « %I. (B8) 

Based on PDF theory (Reeks 1993), the second term on the right-hand-side of (B7) is 
equivalent to the particle velocity covariance multiplied by —2 t~ 2 , while the third term 
equivalent to the particle velocity covariance multiplied by r~ 2 . 

For Sv u'/Ur], we can simplify the relations for the velocity covariances presented 

in Wang & Stock (1993) to obtain 


(vi(t)vf(t)) _ 2 + StSv(ji/I) 

W 2 _ 2[1 + StSv^/C)} 2 ' 


(B9) 


and 

mrm i 

u’ 2 1 + StSv(rj/P) ■ 


(BIO) 


From these relations, we write the horizontal and vertical particle acceleration variances 
as 


and 


(air _ 

V 

2 

'Sv' 


2 StSv + 3 (j) 

a 2 

. u v _ 


St 


2 (stSv+^y 


(ai) p _ 

V 

2 

'Sv' 


1 

a 2 

. u v. 


1st 


StSv + 1 

L v J 


(BH) 


(B12) 


respectively. 


Appendix C. Kinematic collision kernel for an anisotropic particle 
phase 

In this section, we extend the derivation of Sundaram & Collins (1997) to determine the 
kinematic collision kernel formulation for an anisotropic particle phase. From Sundaram 
& Collins (1997), we write the average number of collisions per unit volume over a time 
period r as 


Z c{T) = r ^rf [ ®(r,w,T)p(r,w)drdw, 
" J w J r 


(Cl) 


where n is the particle-number density, $(r,tc,r) is the collision operator, and p(r,w) 
is the joint PDF of the particle-pair separation and relative velocity. At this stage, it is 
simpler to leave the PDF as is rather than split it into a conditional average and the 
RDF as was done in Sundaram & Collins (1997). The collision kernel is given as 


K(d) 


— lim 

n t—> o 


Z c (t) 


2 d 
n 2 dr 


Z c { 0), 


(C2) 
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or equivalently, 

Kid) = [ [ — r ™ H(t*)6[d — r(0)]p(r, w)drdw 

J w Jr HO) 

= [ —H(t*)6[d - r(0)]-^— • f wp(r,w)dwdr (C3) 

Jr nO) J w 

= J —H(t*)S[d - r(0)]^ • g(w p {t)) r dr. 

In the above equations, t* denotes the time required to reach the minimum particle 
separation, H is the Heaviside function, 6 is the Dirac delta function, and g is the 
particle-pair PDF. 

From the definition of t* in Sundaram & Collins (1997), we require r ■ w f 0 for 
H(t*) ^ 0. Since 

r • (w p (t,)) r = ( r p (t ) • w p {t)) r , (C4) 

then in order to satisfy H(t*) f 0, we replace 



• e( wP (t))r, 

(C5) 

with 

-e( wP ~(t))r, 

(C6) 

where 



0 



g(w p ~{t)) r = — J wp{r,w)dw. 

(C7) 


— OO 


Furthermore, since 

^ • Q(w p ~(t)) r = -^e<r p (i) ■ w p ~(t)) r = ^^(r), (C8) 

and we are considering the limit r —> 0 (in which r(0) — > r), we have 


K{d) = f d[d — r]gS^^(r)dr. 
J r 


(C9) 


Since we are interested in collisions on a sphere, we write (C9) in spherical coordinates 
by replacing r with (r, 9, <j>) and dr with r 2 sin 9drd0d(f>. This gives us 

p2,7T pTT 


n Z7T /»7T 

K{d) = J J J S[d — rJpS'Oy (r, 9, cj))r 2 sin9d9d<f>, 


= d 


rJ 0 .70 

/• 27T /* 7T 

2 / / 


g(d, 9,4>)S\ (■ d , 9 , f) sin 9d9d<f>. 


J o Jo 

For an isotropic particle phase, we recover the expected result, 

K(d) = And 2 g(d)S p _ ll (d), 

where 


gW = hJ 


2ir /»7r 


g(d, 9 , (j>) sin 6d9d</>, 


and 


o Jo 

27r />7T 


-1 /* Z7T /»7T 

5^|| (d) = — J J 5^||(d, d, </>) sm9d9d(f>. 


(CIO) 

(Cll) 

(C12) 

(C13) 
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For an anisotropic particle phase, however, (C13) no longer holds. Note that S^(d) is 
an average over particle pairs on the entire surface of the sphere, while S^n (d, 9, cj>) is 
an average over particle pairs on a differential element of the sphere. When the particle 
phase is anisotropic, we must compute S^||(d) by weighting each of the averages on the 
differential elements 5^,, (d, 6, (j>) by the number of particle pairs contributing to that 
average. In this case, we have 

S P _ l{ (d) = Ir C £ S -\\ {d,9 ^ ) 9{d g(d)^ Sin9d9d(C14) 

which reduces to (C13) for an isotropic particle phase. From (C14), we can therefore 
show that 

/* 27T /»7T 

4nd 2 g{d)S p _ [l (d) =d 2 J J g(d,9,cj))S p _ ll (d,6,<j>)sm6ded(i), (C15) 

which is precisely the result we derived in (CIO) for the anisotropic collision kernel. 
We have therefore demonstrated mathematically that (Cll) holds for both isotropic 
and anisotropic particle phases. This explains the empirical observations in Ayala et al. 
(2008), from which the authors argue that (Cll) holds even when the particle distribution 
is anisotropic. 
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